NASA TECHNICAL NOTE 




L0&N COPY: R£ = 
AFWL (W°== 
Kl^tlAND AFI «r — 


A FAST METHOD OF 
ORBIT COMPUTATION 

by 

Karl Stumpff 

Goddard Space Flight Center 
Greenbelt, Md. 

and 

. :4a 4 / 

E. H. Weiss / v r A; 

4 

ZBAf Federal Systems Division r < 4 •; . 

Greenbelt, Md. 




NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


WASHINGTON, D. C. 


APRIL 1968 


TECH LIBRARY KAFB, NM 


tech library kafb, nm 

■■mini 

0l31DSb 

NASA TN D -44 / u 


A FAST METHOD OF 
ORBIT COMPUTATION 


By Karl Stumpff 

Goddard Space Flight Center 
Greenbelt, Md. 

and 

E. H. Weiss 

IBM Federal Systems Division 
Greenbelt, Md. 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


For sole by the Clearinghouse for Federal Scientific and Technical Information 
Springfield, Virginia 22151 — CFSTI price $3.00 




ABSTRACT 


The problem of rapidly computing trajectories of spacecrafts from 
their initial conditions has become very important. Classical methods 
rely almost exclusively on precise integration techniques, but results 
thus obtained are too slow over extended arcs, even on high-speed com- 
puters. Moreover, great accuracy is often unnecessary. Here we pre- 
sent a new method of computing approximate ephemer ides of a small body 
(minor planet or artificial satellite). This method is 10 to 15 times 
faster than the well-known methods of Encke or Cowell. The errors are 
small (e.g., of the order of one part in a thousand) and the results con- 
verge to theN-body point-mass solution for small time steps. It is also 
possible to account for non-point-mass effects; this, however, has not 
yet been implemented. 
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A FAST METHOD OF ORBIT COMPUTATION 


by 

K. Stumpff* and E. H. Weiss** 
Goddard Space Flight Center 


INTRODUCTION 

The new method of computing perturbations described below has been developed in response 
to a need for a means of quickly approximating a spacecraft ephemeris. In particular, a quick 
computation of the orbit of Explorer 33 (described in Example 3 of ’’Results") was required. The 
results need not be exact; small errors, of the order of one part in a thousand, are permitted. 

The special perturbation method yields: 

a. The fast exact solution, or an even faster approximation, to the N-body 
point-mass problem. 

b. Fast approximations to many non-point- mass problems. 

The exact solution to the non-point- mass problem can be obtained by numerical integration. This 
solution, however, has not yet been implemented. 

Reference 1 describes a forerunner of the present technique, used as early as 1942 to com- 
pute heliocentric orbits of minor planets under the influence of Jupiter and other major planets. 
The proof in this report is shorter, more direct, and more convenient for modern applications 
than its earlier counterpart. Both methods are variants of the well known Encke special pertur- 
bation technique. Encke’ s perturbations are defined as the deviations of the planet’s coordinates 
from those of an osculating Kepler ellipse and are of the order h 2 , where h = t - t 0 is the inter- 
mediate time beyond t 0 , the epoch of osculation. The present technique (as in Reference 1) com- 
bines several Keplerian orbits to form an intermediate orbit that includes essential parts of 
Encke’ s perturbations. The deviations of the actual from the intermediate orbit are termed "rest 
perturbations." They are of the order h 4 and therefore very small for small h. Encke’ s per- 
turbations and Stumpff’ s rest perturbations cannot be solved in closed form; they are solved by 
classical numerical integration. 


*Prof. Emeritus, Univ. of Goettingen and Senior Post-Doctoral Resident Research Associate, U. S. National Academy of Science- 
National Research Council. 

**Advisory Mathematician IBM Federal Systems Division, P. O. Box 67, Greenbelt, Maryland. 
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NOTATION AND EQUATIONS 


The results about to be derived are valid for N point-mass bodies. For simplicity the formu- 
las are given for the 4-body case (e.g., earth, moon, sun and spacecraft), but the proof is readily 
extended to N bodies. The formulas exhibit remarkable symmetry. 

Let each of the subscripts i, j, k, and 1 assume the value 0, 1, 2, or 3 with the proviso that dif- 
ferent subscripts be distinct. Denote the four bodies by Q. and their masses by m. . The vectors 
from Q i to Q. are q.., and their magnitudes are r.. = |q ij |. Without loss of generality, we can 
use canonical units; then the Gaussian constant equals unity. Assume that the bodies act as point 
masses. Then the vectors q.. satisfy the differential equations 


-0"i + m j) Q i j r i j ” 3 ~ m k( q ik r ik 3 + 




+ q ii r !j 


')• 


(1) 


The 12 combinations i, j from 0, 1, 2, 3, contain only three linearly independent vectors q £j (for 
instance q 10 , q 12 , q 13 ), as there exist six identities q ^ = -q.. , and three independent equations 
of the form 


q i j + q jk + 


o. 


Here and elsewhere, q .. (0)= q.. (t 0 ) refers to the time t 0 ; q £j = q.. (t) refers to the time 
t = t 0 + h . A special solution of Equation 1 is determined by the initial values 

q i j (0) andq ij (0). (2) 

Use square brackets to denote Kepler ian orbits that osculate q^ at t = t 0# Then 


[qij(O)] = Qij (0); [q ±j (0)] = q i j ( 0 ) 


(3) 


are the conditions of osculation. 

The osculating Keplerian orbits satisfy the differential equations 

[qi.] = _( m . +m,) [q^ ] [r^]' 3 , (4) 

where 

Cr , , ] = | [q i j ] I • 


Equation 4 can be applied to any combination of two different subscripts i, j, k, and 1. Now 
define the vectors s.. (=-s..) by 


[q u ] [r j j ] -3 


(5) 
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Substituting Equation 5 in Equation 4 and rearranging, one obtains 


i r ii ~ 3 = — - — ( 6 ) 

1J1J m.+m. J J 

i j 

Using Equation 6 to eliminate all expressions of this form from the right side of Equation 1 gives: 


q = [q ] + 

1 J 1J m . + m. 


tq ik ] + 


*k T “‘J 


+ 


m . + m . 


[<3 U ] 


Cqij] + ( m i + m j ) s ij + m k < s ik + s kj) + m i < s x i + s ij )• 


This can be written 


q u = [q i tS.j + R ij 


(7a) 


where 


S. . ^ 


(q.J 


m . m . 

(q k J + [q u ( + i 

m . x m . 11 m 4- m . 1 J 


R. 


ij = ( m i + m j ) s ij + m k < s ik + s kj) + m i ( s ii + s ij )• 


(7b) 

(7c) 


The first and second integrals of [q^ J and S. . can be obtained from the well known properties of 
Keplerian conic sections. See Reference 2, or Chapter V of Reference 3. The integration of L, 
however, can be effected only by numerical methods. The first and second integrals of Equation 
7a are 


= ld >i j ] +s tj +Rij + a ij- 

q i j = c«iij 3 +s ij +R ii + a ii h + b ii 


(8) 


To determine the constant vectors a., and b . . , evaluate Equation 8 at t 0 , from which, by 
Equation 3, 


S ij (0)+R ij (0) + a. j = 0 , 
S ij (0) + R ij (0) + b ij = 0 . 


Without loss of generality, (0) and R^ (0) can be equated to zero, since non-zero values can be 
absorbed by a., and b.. . Thus 


= - s ij(°)- b ii 


(9) 
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Substituting Equation 9 in Equation 8 gives 


where 


= [ ^ii ] + Pij + Kij 
q ij = t q ij J + p ij + R ij * 

Pi. =S U -Si, (0). Pi, =s ij -h Si,(0) -Sj, (0). 


( 10 ) 


( 11 ) 


P Aj is termed the approximate perturbation, and r. . the rest perturbation. Equation 5 shows 
that s.. (0) = o, whence R^ (0) = 0 in view of Equation 7c. Similarly r (0) = 0, as can be seen by 
differentiation, and hence R.. is of the order h 4 . For sufficiently small values of h, R.. and r.. 
in Equation 10 can be ignored, if desired. The resulting approximations are 


( 12 ) 



Exact results are given by the formulas 



(13) 


where q and q * j are obtained by Equation 12, and R^ and R^ are the first and second inte- 
grals of R. j , defined by Equation 7c. 


USAGE 

There are three different ways to compute (accurately or approximately) the ephemerides of 
celestial bodies by the set of formulas just given. For computational simplicity, one of the four 
bodies is placed at the center of the coordinate system. This body is denoted by Q c and is called 
the central body. Let r be a dummy subscript that can assume the values 0, 1, 2, or 3 as long as 
c / r. Then the initial conditions are given by the known values q cr (0) and q cr (0) , and our prob- 
lem is to find q cr and q cr . By vector addition any other vector between the four bodies can be 
determined since 


where 


^ i j ^ i c + q cj ’ 
^ i c ~ ** ^ci 
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Method A 


Method A, described in Reference 1, is a variant of Encke T s method of computing special per- 
turbations. This method is especially effective when the numerical integration uses difference 
tables; this mode of integration was especially popular before the advent of high-speed computers. 

Compute the six osculating Keplerian orbits [q . .] and their velocity vectors [q..] at t = t 0 +nh 
( n= l 9 2 , . . , N) ? where h is a conveniently chosen constant step length. Then compute S cr (t) and 
S cr (t), the two-body solutions of Equation 7b satisfying the initial conditions given by Equation 3. 
Compute P cr (t) and P cr (t) by Equation 11, and obtain the coordinates of the intermediary orbits 
q* r (t) and q* r (t) by Equation 12. The rest perturbation R cr (t) can be obtained by integrating 
Equation 7c twice, using the classical method of numerical integration by differences, with R(0) = 

R (0) =0 as initial conditions. The start of the integrating scheme of differences can be achieved 
without iteration: Calculate s .. (t) from 


S*. = [q..] [r..]“ 3 — q ? . r*. ~ 3 
i J M i j i j 4 ij i j 


(14) 


instead of Equation 5 for small nh # Since s.. - s*. is of order h 4 it follows from Equation 7c 
that the error of R i} is also of order h 4 . Therefore the error in R.. is of order h 6 and can be 
neglected for small intermediate times. This is a remarkable advantage over Encke’s method, 
where iterations at the start cannot be avoided. 

Method A is effective if (1) none of the bodies closely approach any other body, (2) the mag- 
nitudes of R remain substantially smaller than the magnitudes of the complete perturbations, P+R. 
Thus the calculations are restricted to relatively small time intervals. 


Method B 

Let h 0 = tj - t 0 be a small step length, for which R(tj ) is very small. Compute 
^cr (t i ) ~ ^cr ( t i ) 9 ^cr ( 1 1 ) ~ ^ c r (* 1 ) by Method A, neglecting the rest terms. Use these approx- 
imations of q c t x ) and q (t t ) as new initial conditions, and compute q(t 2 ), q(t 2 ) at t 2 = t 1 + hj, 
again neglecting the rest terms, etc. The results will be accurate as long as the neglected rest 
terms do not accumulate beyond a specified tolerance. Thus Method B, which is simple and effi- 
cient, can be used to approximate orbits when great accuracy is not needed. 


Method C 

Use Method B to compute q* r (t 0 + 4h 0 ) as well as q* r (t 0 + nh 0 ) for n = 1, 2, 3, 4. Then, 
using s*. defined by Equation 14 instead of s.. , compute R cr (t 0 + nh 0 ) by Equation7c, and obtain 
R (t o + 4h 0 ) and R(t 0 + 4h 0 ) by integration. For example, Stirling’s five-point formula with 
R (0)= 0, and with the initial conditions R(t 0 ) = R (t 0 )= 0, yields 
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(15) 


R(t 0 + 4h 0 ) = ^| [64R(t 0 +h Q ) + 24R (to + 2h 0 ) + 64R(t 0 + 3h 0 ) + 14R(t 0 + 4h 0 )] 




h 0 2 


R(t 0 + 4h Q ) = [l92R(t 0 +h 0 ) +48R(t 0 + 2h 0 ) +64R(t 0 +3h 0 )]. 


Then use Equation 13 to obtain q (t t ) = q* (t t ) + R(t , q(t = q*(t + R(t at t x = 1 0 + 4h 0 . 
With these vectors as new initial conditions proceed tot 2 = t 1 + 4h 1 , t 3 =t 1 + 4h 2 , etc. 

This method delivers accurate ephemerides for an extended range of time as long as the step 
lengths h n are sufficiently small, though they may be noticeably larger than those used in Method 
B. If the h n surpass a certain limit, errors will occur because: 

1. The supposition s « s* is no longer valid. This may be cancelled by iteration, if desired. 

2. Integration Equations 15 use interpolating polynomials of the fourth degree. The formulas 
lose accuracy for large h, but this error can be avoided by the use of integration formu- 
las of a higher degree. 

The relative smallness of the rest perturbations permits the use of rather large time steps. 
Proper time- step control significantly reduces the computing time; considerable effort was spent 
in selecting the "best” time-step. The four time-step criteria given below were specialized to the 
following bodies and units: 

1. Q 0 : spacecraft 
Q x : earth 

Q 2 : moon 
Q 3 : sun 

2. Length is measured in earth radii, time in canonical units of 806.813645 seconds, and the 
mass of the earth is taken as unity. The four time-step criteria are: 

a. h is constant. This is best avoided unless the orbit is nearly circular. 

b. h = (K/Q) 1/4 ? where K is an input parameter, and Q (a theoretical overestimate of the 
rest perturbations) is given by 



3. h = (A + W)/B , where W = Min (r 10 ? Cr 20 ) ? and A, B, C are input parameters. 

4. h new = (l/2)( |d/V 3 4 R 10 1 +i) h old , where D is an input parameter, and 

V 4 R 10 =R 10 (t 0 +4h) -4R 10 (t 0 + 3h ) + 6R 10 (t 0 + 2h) -4R 10 (t 0 + h ). 

Large changes in step size are avoided by the condition 0.5 h old < h new < 3.5 h old . 
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The following section gives numerical examples. In Tables 3, 4, and 5, the column headed 
"time-step criterion" contains either: the fixed value of h; or K; or A, B, C; or D. In general, best 
results are obtained with the A, B, C, or D criteria. 

Test results were computed not only with large step sizes but also with very small ones. It 
was always possible to converge to results that are invariant under further reduction of the time- 
step. These converged results may rightly be considered the accurate ephemerides of the 4-body 
point-mass problem. 

RESULTS 

The following five examples illustrate the results. 

Example 1 

Method A was first used by K. Stumpff (Reference 1) in 1942 to compute special perturbations 
of the minor planet (931) Whittemora by Jupiter. Some numerical results are taken, abbreviated, 
from Reference 1, to illustrate the difference between this method and the classical method of 
Encke (see Reference 4, p. 378ff). The formulas of Section 1 may be used with masses 

m 0 = 0 (minor planet), 

inj = 1 (sun), 

m 2 = 1/1047.35 (Jupiter), 

m 3 = 0 (as no other perturbing planet has been taken in account). 

The initial epoch is t 0 = 1920, April 29.0; and h, the constant step of integration, equals 
40 days = 40k canonical units (k - 0.0172021). 

Table 1 gives the heliocentric equatorial x-coordinates of the minor planet for t = t 0 + nh 
( n =0, 1, 2, 8). It lists 

P = approximate perturbation (Equation 11), 

R = rest perturbation obtained by numerical integration of Equation 7c, 
p+R = complete perturbation, 

= complete perturbation derived by Encke's method. 

Table 1 illustrates: 

1. The very slow increase of R for small intermediate times, compared with the 
rapid increase of P and cr, 

2. cr = p + r, except for small deviations due to rounding. 
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Table 1 


Example 2 


Perturbations of the x-coordinate of (931) 
Whittemora for t 0 + 40 n Days. 

n P R P+R cr 

0 0.0 0.0 0.0 0 

1 8.6 0.0 8.6 8 

2 38.5 0.4 38.9 39 

3 95.5 3.9 99.4 99 

4 182.9 17.9 200.8 201 

5 301.4 55.0 356.4 356 

6 448.8 133.4 582.2 581 

7 620.8 276.6 897.4 896 

8 807.8 513.6 1321.4 1320 

NOTES: 1. t 0 = 1920 April 29-0. 

2. Results expressed in 10 A.U. 

Table 2 

Coordinates of (931) Whittemora at 
t 0 + 1600 Days, in A.U. 

x y z 

0.6061220 2.299933 -0.2915877 

0.6061201 2.299943 -0.2915872 

0.6061180 2.299950 -0.2915866 

0.6061054 2.299979 -0.2915840 

0.6060915 2.300004 -0.2915813 

0.6060304 2.300097 -0.2915702 

0.6061222 2.299931 -0.2915879 

NOTES: 1. t Q = 1920 April 29-0. 

2. Results expressed in A.U* 


h in days 

20 

40 

50 

80 

100 

160 

Exact Value 


Method B is used to compute the orbit of 
the same minor planet, (931) Whittemora. 

The computations extend to 1600 days (ap- 
proximately 80 percent of one revolution) 
beyond the epoch 1920, April 29.0. Table 2 
lists the terminal heliocentric rectangular 
coordinates of the planet in A.U. for several 
constant values of h , as well as the true values 
for comparison. The true values are obtained 
almost perfectly when h = 20 days, and the 
error barely exceeds 10 -4 A.U., when the 
computation is performed in 10 steps of 160 
days each. 

The remaining examples involve the ap- 
plications of Methods B and C to a spacecraft 
in the gravitational field of earth, moon and 
sun. The orbits are highly eccentric, closely 
approaching the earth or the moon. Hence, 
variable step lengths should be chosen. 

Example 3 

Compute the orbit of Explorer 33, 
launched July 1, 1966. Integrate for 180 days 
beyond t 0 , the epoch of computation, where 
t 0 = 1966, July 31. The spacecraft describes 
over 12 highly eccentric trajectories around 
the earth and several times closely ap- 
proaches the surfaces of the earth and the 
moon (to within approximately 3 and 5 earth 
radii, respectively). The earth- spacecraft 


separation is the criterion that measures the effectiveness of the method. The exact 4-body point 
mass separation is 138801.68 km at t 0 + 170 days and 437108.24 km at t Q + 180 days. Table 3 
gives the deviations of the earth- spacecraft separations at these times, for several runs involving 
different time- step controls. The deviation always attains its maximum near t Q + 170 days. 

Table 3 also gives the machine execution time and the number of steps used in the computation. 

The last line gives the equivalent information for the JPL-Holdridge program (Reference 5) which 
accounts for planetary and non-point-mass perturbations. The agreement, it will be noted, is quite 
satisfactory. 
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Table 3 


Explorer 33: Deviations in Earth Spacecraft Separation, 


Deviation of earth- 
spacecraft distance 

Deviation of earth- 
spacecraft distance 

Time 

No. of com- 

Method 

Time -step 

at t 0 +170 days 

at t 0 + 180 days 

(seconds) 

puting steps 

criterion 

(km) 

2495 

(km) 

915 

50 

Not Available 

B 

K = 10* 5 

766 

279 

88 

Not Available 

B 

K = 10' 6 

245 

90 

157 

Not Available 

B 

K = 10 -7 

41 

11 

421 

2440 

C 

4h = 3 hrs 

8146 

3052 

38 

218 

C 

A = 0 

940 

286 

53 

310 

C 

B = 1.5 
C = 1.5 

A = 0 

722 

220 

61 

356 

c 

B = 2.5 
C = 3.5 

A = 0 

616 

186 

70 

408 

c 

B = 2.5 
C = 1.5 

A = 0 

86 

27 

86 

502 

c 

B = 2.5 
C = 1 

A = 0 

42 

13 

122 

711 

c 

B = 4 
C = 2.5 

A = 0 

1 

1 

190 

1104 

c 

B ~ 4 
C = 1 

A = 0 

0 

0 

363 

2108 

c 

B = 8 
C - 1.5 

A - 5 

930 

436 

about 15 min 

Not Available 

JPL-Holdridge* 

B = 16 
C = 1 


* - Accounts for planetary and non -point mass effects. 


The inverse relation between speed and accuracy can be seen best from Figure 1. It plots 
maximum deviation of the earth- space craft distance versus machine time for Explorer 33, com- 
puted 180 days beyond the epoch. Note that Method C is superior to Method B in this example and 
probably in others. 
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NOTES 



COMPUTING TIME (sec) 


Figure 1— Maximum deviation in earth-spacecraft 
separation for Explorer 33. 


1. All machine times refer to the IBM 7094 
Model 1 . 

2. The number of time -steps is available for 
Method C only. A time- step extends 
from t.tot i + 1 =t.+4h i . 

Example 4 

Consider a typical earth- moon trajectory 
integrated for 2.5 days beyond the epoch. The 
trajectory starts less than 450 km from the 
surface of the earth, i.e., at third-stage cut- 
off. The exact earth-spacecraft separation 
for the point- mass earth, moon, sun problem 
equals 362,141.35 km at t 0 + 2.5 days. Table 
4 lists the effects of different time -steps 
and includes the results of the JPL-Hold- 
ridge program, which accounts for planetary 
perturbations. The results agree very well 
with the point-mass JPL-Holdridge program. 


Example 5 

Consider a typical lunar orbit with the following initial characteristics: aposelenium = 8634 
km, periselenium = 3443 km, period = 702 min, e = 0.429. The integration extends for 180 days 
beyond the epoch July 4, 1966. Table 5 gives the terminal discrepancy in moon- spacecraft separa- 
tion for different time-steps. The exact value is 8562.37 km. There is close agreement between 
the simple 4-body solution and the JPL program, even though the latter accounts for harmonics 
and planetary perturbations. 


COMPUTER PROGRAMS 


Four double precision computer programs have been written in FORTRAN IV. They are en- 
titled STUMPFF1, STUMPFF2, STUMPFF3, and STUMPFF4, and are reproduced in Appendix A. 
This section briefly describes points of interest to users of the programs. 

The programs generate their own ephemerides; this facilitates programming and saves 
memory locations. 
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Table 4 

Typical Earth- Moon Trajectory: Deviations in Earth- Spacecraft Separation at t 0 + 2.5 Days. 


Deviation of earth- 
spacecraft distance 
(km) 

Time 

(seconds) 

No. of 
steps 

Method 

Time-step 

criterion 

Comments 

0.73 

1.44 

7 

C 

D = 5 x 10- 6 


0.37 

1.63 

8 

C 

D = 2 x 10' 6 


0.06 

1.99 

10 

c 

D = 9 X 10' 7 


0.03 

2.18 

ii 

c 

D = 7 X 10" 7 


0 

3.26 

17 

c 

D = 2 X 10" 7 


0.51 

about one 
minute 


JPL-Holdridge 


Planetary perturbations 
included. Point-mass 
bodies assumed. 

2791.98 

about one 
minute 


JPL-Holdridge 


Planetary perturbations 
included. Earth and moon 
harmonics included. 


STUMPFF1 


The program is set up to compute Example 2 by Method B. The central body is the sun; lo- 
cations XO, YO, ZO, XDO, YDO, and ZDO contain the initial position and velocity values of (931) Whitte- 
mora; XI 0, Y10, Z10, XD10, YD10, and ZD10 contain the initial values for Jupiter. The unit of 
mass is the mass of the sun; the mass of Jupiter = 1/1047.35; the mass of the minor planet equals 
zero. The unit of length is the A.U. and the unit of time 58 d . 13244. The step size, in days, is in 
location DIFF. The program prints the number of days beyond the epoch, and the coordinates of 
the minor planet and Jupiter. The program stops when TAU > TAUMAX, where TAU = 0.0172021 
X DIFF. 

STUMPFF2 


The program is set up to compute Example 3 by Method B. Q 0 is the spacecraft; Q x is the 
earth, which is the central body; Q 2 is the moon; and Q 3 is the sun. Locations Y10(I), Y12(I), and 
Y13(I), (I = 1, 2, 3, 4, 5, 6), contain the initial values of q 10 (0), q 10 (0), q 12 (0), q 12 (0), q 13 (0), q 13 (0) 
in canonical units. The unit of length is the mean earth radius of 6378.165 km, the unit of mass is 
the mass of the earth, and the unit of time equals 806.813645 seconds. The program prints the 
initial conditions. Then it prints four lines every N th day, namely: 

(Line 1) q 10 , q 10 , number of days since epoch, 

(Line 2) q l2 , q 12 , 

(Line 3) q 13 , q l3 , 

(Line 4) q, r 0 , r 20 , number of days since epoch. 
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Table 5 


Typical Lunar Orbit: Deviations in Moon- Spacecraft Separation at t 0 + 180 Days. 


Deviation of moon- 
spacecraft distance 
(km) 

Time 

(minutes) 

No. of 
steps 

Method 

Time -step 
criterion 

Comments 

1068 

6.9 

2416 

c 

A = 2 
B = 1.5 
C = 1 


43 

13.8 

4814 

c 

A = 2 
B = 3 
C = 1 


51 

17.5 

6104 

c 

A = 5 
B = 5 
C = -1 


14 

23.0 

8020 

c 

A = 2 
B = 5 
C = 1 


22 

25.2 

8762 

c 

A = 6 
! B = 9 
C = -1 


0 

87.5 

30459 

c 

A = 2 
B = 19 
C = 1 


28 

29.8 

9230 

c 

D = 3 X 10' 7 


16 

35.8 

11108 

c 

D = 2 X 10' 7 


0 

68.1 

21108 

c 

D = 6 X 10' 8 


28 

about 120 


JPL-Holdridge 


Includes harmonics and 
planetary effects. 
Encke mode used. 

22 

about 175 


JPL-Holdridge 


Includes harmonics and 
planetary effects. 
Cowell mode used. 


The results are printed in km and km/sec. At present, the program prints every tenth day. This 
can be changed by altering "10.D00" in the two consecutive instructions: 

IJI = IDINT (TIMED/ 1 0.D00), 

IJ2 = IDINT (TIMEDN/ 1 0.D00). 

The program stops TIMEMX days beyond the epoch. 


STUMPFF3 

The program computes Example 5 by Method C, using the A, B, C time-step criterion. The 
bodies are numbered as in STUMPFF2 and the units of length, mass, and time are defined as in 
STUMPFF2. The first printout gives the initial conditions followed by A, B, C, and m im Subse- 
quent printouts are four lines each: 
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(Line 1) q 10 , q l0 , number of days since epoch 
(Line 2) q 12 , q 12 
(Line 3) q 13 , q 13 

(Line 4) Contains five words. The first is immaterial; the others are r 10 , r 20 , number of 
days since epoch, and number of computing steps. 

There are eight input cards per case, and two or more cases may be stacked. The fields of 
the input cards end in columns 16, 32, 48, and 64, and contain the following floating point input: 

(Card 1) q i0 (0) 

(Card 2) q 10 (0) 

(Card 3) q ia (0) 

(Card 4) q 12 (0) 

(Card 5) q 13 (0) 

(Card 6) q 13 (0) 

Card 7 has four fields that specify A, B, ONEDAY, and TIME MX. The only field of Card 8 speci- 
fies C. The time-step criterion uses A, B, and C; ONEDAY governs the frequency of the printing; 
and a case is terminated after TIME MX days of computation. If column 72 of Card 7 equals 1, a 
new case will be processed after the present case; the last case must contain a blank in column 72 
of Card 7. 

STUMPFF4 


The program computes Example 5 by Method C, using the D time -step criterion. Everything 
is as in STUMPFF3 except: 

1. Input Card 8 does not exist. Card 7 contains D, h init , ONEDAY, TIMEMX. D governs the 
time step; h init gives the initial value of h; and ONEDAY and TIMEMX are as in 
STUMPFF3. 

2. The line printed after the initial conditions contains D, h min , one immaterial word, and m 1 . 


Goddard Space Flight Center 
National Aeronautics and Space Administration 
Greenbelt, Maryland, October 2, 1967 


311-02-01-01-51 
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Appendix 

Symbolic Listing of FORTRAN IV Programs 
STUMPFF1, STUMPFF2, STUMPFF3, and STUMPFF4 

$ JOB 0510 S TUMP F F 1 BY E.H. V%i & I S S RIGGS BLCG X 7319 

SPAUSE 

$E XECUT E I B JOB 

$ I B JOB GO, LOGIC 

C STU W PF F 1 BY E.H. WEISS RIGGS BLCG X 7319 

$ I B FTC VAIN L I S T , REF , ChCK 

DOUBLE PRECISION 

DGK,OIFF, TAUMAX, TAU, T AU1 , XO, Y 0 , Z G ,XCu, YDC, ZOO, 
DX10,Y1C,Z10,X010, YDI0,Z010,CKAP0 ,QKAP 1 , GKAP2 , 
DTE1,TE2,TE3,TE4,TE5,TE6,TE7,TE3,TE9,TE1C,TE11,TE12, 

DSX, ^Y, SZ ,SCX f SDY, SUZ » S IGX » SI GY, S IGZ 
QK= 0 . C 1 72 02 1 DOO 
0IFF=100.CC0 
WRITfc (3,30) 

30 FORMAT ( IHl) 

T AU W AX = 2 3 . DCO 

TAU= O.DCC 

TAU1=QK*CIFF 

X0=-3. 244464000 

Y0=.28826CDC0 

Z0=. 5726130C0 

XDO*-. 164671D00 

YDO = - . 5064 19000 

ZD0=.069948DC0 

X10=-4.09619DC0 

Y10=3. 11893DG0 

Z10=1.439800C0 

XD 1C = - . 19726DC0 

YD10=-. 2013b 000 

ZD10=-. 08162000 

wK A °0 = l.CDCO 

CK A P 1 = 1 .CDC0/1047. 35UCG 

Q K A P 2 = 1 .CDC0 + GKAP1 

TE1=DSGRT { GKAP2 )/ { G K * 4 C .000 ) 

XDIC=XD1C*TEI 
YD 1 P=YD1C*TE1 
ZD10=ZC1C*TE1 

31 T AU = TAU+ TAU1 

CALL SUSM TAUlfXOfYO,ZOf XDC, YDO , ZOO , OKAPC , 

1 SX,SY,SZ,SCX,SDY,Sl;Z,SIGX,SIGY, SIGZ) 

TE 1 = XC+SX 
TE2 = YO+SY 
TE 3= ZC + SZ 
TE4 = XCO+SCX 
T E 5 = YCO + SDY 
T E6 = ZDO+SCZ 
TE 7 = X0- X 1C 
TE8 = Y0-Y 10 
TE9=Z0-Z 1C 
TE 1 0= XDO-XD 10 
TE 1 1=YC0-Y010 
TE 1 ?=ZD0— ZC 10 

CALL SUB1( TAU1,TE7, TE8, TE9,TE10, TE11, TF 12,QKAP1, 

1 SXf SY«SZ« SOX vSDYfSUZvSr GXvSIGYv SIGZ) 
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TE1=TE1+SIGX 
TE2 = T E2+ $ I G Y 
TE3 = TE3+SIGZ 
TE4=TE4+ SDX 
T E5 = TE5+ SDY 
TE6=TE6+ SDZ 

CALL SUBi(TAUlfX10*Y10t ZlOfXDIO, YC10*Z010,QKAP2, 

1 SX,SY,SZ, SDX, SDY, SUZ,SIGX ,SIGY, SIGZ) 

XO= TE1+ CKAP1*SIGX 
YO= TE2+ CKAP 1* S I GY 
ZO=TE3+QKAPl*SIGZ 
XDG = T E4 + G K AP 1 * SDX 
YDO = TE 5+ CKAPi*SDY 
ZDO = TE6 + CKAP 1 *SOZ 
X 1 0= X10+SX 
Y 1 0 = Y1C + SY 
Z10= Z 10+SZ 
XD 1 0 = XD 1 0+ SDX 
YD 1 0 = YD 1 C+ SDY 
ZD10=ZD10+S0Z 

WRITE ( 3,35 ) TAU,XO, YO,ZO ,XCO, YDO, ZOO 

35 FQRVAT ( 1X6H TAU =,D1S.7,13H S j v A L L BCD Y = , 6 ( D15 • 7 , IX ) ) 

WRITE (3,36) X1C, Y1G,Z10, XD 10, YD 10, ZD10 

36 FORMAT ( 1 X34H JUPITER= . 6 < Cl 6 . 7 , IX ) > 

IF(TAU.LE.TAUMAX)GC TG 31 

RETURN 

END 

$ I BF TC SUR1 LIST, REF ,CtCK 

SUBROUTINE SUB1 (TAU,XO,YO,ZC,XDC, YDO, ZOO, CKAP, 

1 SX,SY,SZ,SDX,SCY,SDZ,SIGX,SIGY,SIGZ) 

DOUPLE PRECISION 

D TAL',XO, YC,ZO,XDO, YOC , ZCO , CK AP , 

D SX,SY,SZ,SDX,SDY,SDZ,SIGX,SIGY, SIGZ, 

D RO,CMUO,CKSIO,SIGO,ETAO,OWEO,EPSO,ZETAO,RHOO,CHIO, 

DTEMP1 , TENP2, 

DZ, DELTA, Cl, C2,C3, 

D FM,G,FD,GDP 

RO-CSQRT ( X0**2+Y0«*2+Z0**2) 

QMUO= GKAP/R0**3 
QKSrO= QML 0 * T AU* * 2 

S I G0= ( X C*XDO+ YO * YDO + ZO*ZCO ) /RC**2 
ETAG= S I GC* T A U 

OME 0= ( XD0**2+YD0**2+ZD0**2 ) /R0**2 
E PSC- OMEC- QMUO 
ZETA0=EPSC*TAU**2 
RHO 0 = CMUO- EPSG 
CH I 0= RHCC*TAU**2 

CALL SUB2(ETA0*ZETA0,CHI0t Z»CELTA,C1. C2,C3) 

FM= C 2 *QKS I 0* Z* # 2 
G=TAU*(1.D00-(C3*QKSI0*Z**3) ) 

FD=-(Cl*QKSIO*Z)/(DtLTA*TAU) 

GDM= C2*QKSI0*Z**2/UELTA 

SX= G* X DO- F M * X 0 

SY= G* Y DC- FM* YO 

SZ= G * Z D C- F^*ZO 

SDX = FD* XO- GDM*XDO 

SDY= FC* YO- GDK*YOG 

SDZ = FD* ZO- GDM*ZDJ 

TEMP1=-QKSI0*Z**2 

TEM P2 = C 3* Z * T AU 

S I GX= TENP1*(C2*X0+TEMP2*XC0) 

S I G Y= TE. V P1*{C2*Y0+TEMP2*YC0) 

SIGZ=TEMP1*(C2*Z0+TEPP2*ZD0) 

RETURN 

END 
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SIBFTC SUB2 L I ST , REF , DEC K 

SUBROUTINE SUB2 [ E T A , l E T A , CH I , Z , CE L T A f C 1 , C2 # C 3 ) 

DOUBLE PRECISION 

D QL A 9 Z, Cl, C 2 , C 3 , ETA, ZETA, CHI, 

DDELTA,TOLl ,CH 
ITER=C 
Z = 1 * ODOO 
TOL = 1 « D-08 
201 ITER=ITER+1 
QL A=CH I * Z * Z 

C2=.5DOO*( 1 .000- ( GLA/12.D0C)*( 1. COO- ( QL A/ 30 . DC 0 ) * 

1 { 1 .DCG-(QLA/56.000) »< l.DCO-l QLA/90.C00) >)) ) 

C3= ( l.DC0/6.D00)*( 1 . 000- < QL A/20 . DOC )* ( 1 . 000- ( QL A/42 . DOO ) * 
1 <1.D00-(QLA/72.D00)*(1.DCC-(QLA/110.D00> ) ) ) ) 

C1=1.DC0-(QLA*C3) 

DELT A= 1 . DC0 + C1*ETA*Z + C 2*Z ET A*Z * *2 
QH = C2*ETA*Z**2+C3*ZETA*Z**3+Z-1.DC0 
Z = Z- QH/CELTA 

I F ( DABS ( QH ) .LE.T0L1 ) RETURN 
I F ( ITER.LE. 1 0 ) GO TO 201 
RETURN 
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$ JOB 0510 STUMPFF2 BY E.H. WEISS RIGGS BLDG X7319 

SPAUSE 

^EXECUTE IBJOB 

$ I B JOB GG » LOGIC 

JIBFTC MAIN LIST, REF, NODECK 

DOUBLE PRECISION 
DERRtTIMECN, ITDN, ITD, 

DX10,X12,X13.XTEMP,XKEP,XDEL, 

DM1,M2,M3,M12,M13,M23,M2F12,M2F23,M3F13,M3F23, 
DCML,CMV,CMT, Y10, Y12, Y13,TIMED,TIMEMX, 

DTAUO, TAU,R10,R20,TEMP,Q,R10E2,R2CE2, 

DZ1,Z2,Z3,Z4,Z5,Z6, 

DL ,C1 iC2 ,C3, DELTA, H, 

DRE2,R,KSI,ETA,ZETA,CHI 

COMMON X10(12)»X12(12)»X13(12)*XTEMP(6) ,XKEP(6),XDEL(6), 
CM1,M2,M3,M12,M13,M23, M2F12, M2F23, M3F13.M3F23, 
CCML»CMV,CMT » Y 10 1 6 ) * Y 12 I 6 ) « Y 1 3 ( 6 ) , T I ME D» T l M EMX » 
CTAUD,TAU,R10,R20,TEMP,C,R10E2,R2CE2, 

CZ1,Z2,Z3,ZA,Z5,Z6, 

CL ,C1 ,C2,C3, DELTA, H, 

CRE2.R.KS I.ETA, ZETA.CHI ,DUM{ 10C0) 

C STUMPFF2 BY EH WEISS RIGGS 

Y10tl)=.1835264CD06 
Y 10 (2 >=-.240943 38 006 
Y10 ( 3 ) = — .36452764D03 
Y 10 I 4 ) = . 1C044146D01 
Y10(5)=-.32C81304DC0 
Y10(6)=-.151680OlD00 

Y 1 2 ( 1 ) = .21384734006 

Y12(2)=-.29619053D06 
Y12(3)= - . 16430464006 

Y 12 t 4 ) = . 84600696D00 

Y 1 2 ( 5 ) = . 4 76 26001 DCO 

Y 12 1 6 ) = .17218044DC0 
TIM*=MX=181.0C0 

Y 13 ( 1 ) =— . 93856134D06 
Y13 (2 ) = . 1C949798D09 

Y 13 ( 3 ) = .47486129006 

Y13 (4)=-. 22920520002 
Y13(5)=-.16789843D02 
Y13(6)=-.72817681D01 
CML=6378.165D00 
CMT=806.813645D00 
CMV =CML/CMT 

C INPUT IS IN KILOMETERS ANC SECONDS. 

C MULTIPLY CANONICAL BY CML.CMT OR CMV, TC OBTAIN METRIC . 

C TAUC IS IN DAYS. TAU IS CANONICAL. 

ERR= 1 . D-06 
Q= 1 . D-07 
Ml® 1 • DCO 

M2=1.D00/81.30L50C52DC0 

M3=332951.26580C0 

M12=M1+M2 

M13=M1+M3 

M23=M2+M3 

M2F 12=M2 /M 12 

M2F73=M2/M23 

M3F 1 3=M 3 /M 1 3 

M3F 7 3 = M 3 /M 2 3 

Z 1 = . 9000 
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Z2= • 9DC0 
Z 3= • 9 DCO 
Z4=.9DC0 
Z5=.9DCO 
Z 6= • 9 C C 0 
DO 25 1 = 1,3 
X 10 ( l ) - Y 1C ( I ) /C ML 
X 1 2 ( I ) = Y 1 2 ( l ) / C ML 
X13CI )=Y13( I ) / C ML 
X 10 ( I +3 ) = Y 1 0 ( 1 + 3) /C MV 
X 1 2 ( 1+3) = Y 1 2 ( I + 3 ) / C M V 
XX3 f 1 + 3 ) =Y1 3 ( 1 + 3 ) /C M V 

25 CONTINUE 
TIMED =O.CCO 

WRITE ( 3,26 ) ( Y 10 ( I ) ,1 = 1,6) , TIMED 
hRITE(3,19)(Y12(I ) , I = 1 , 6 ) 

WRITE < 3, 19 ) ( Y13II) , 1 = 1 ,6 ) 

26 FORMAT ( 1H1, 6016*8, D16. 8) 

9 TAU = DSGRT(CSQRT (ERR/O) ) 

TAU=CMIM1( TAU, ICO. DC 0) 
TAUn=TAU-*CMT/8640 0.DC0 
T I M EDN = T IMEC+TAUD 
IJ1=IDINT(TIMED/I0.DC0) 

I J2=IDINT( TIMEDN/ 10.DC0 ) 

I F ( TJl.EG. I J2 ) GO TC 60 
IJ2=IDINT (TIMEDN) 

T I MFDN = I J2 

TAUC=T IMEDN-T IMED 

T AU= TAUD *86400. DCO /CM T 

TIMFD=TIMEDN 

I P R - 1 

61 CONTINUE 

CALL SUBKXlOfMl, Z1 ) 

DO 10 1=1,6 

X10 ( 1+6 ) = XKEP ( I ) 

10 CONTINUE 

CALL SU81(X12,M12,Z2) 

DO 11 1=1,6 

X10U + 6)=X10( 1+6) +M2F12*XDEL( I ) 

X 1 2 ( 1+6) =X K EP ( I ) 

X13 ( 1 + 6 ) =M2F12*XDEL ( I ) 

11 CONTINUE 

CALL SUB 1 ( X 13 , M13 , Z3 ) 

DO 12 1=1,6 

X10 (I +6 ) = X10(I+6)+M3F13*XCEL(I) 
X 12 ( I +6 ) = X12( I+6)+M3F13*XCEL(I) 
X 1 3 ( I +6 ) = X 1 3 ( I+b)+XKEP< I ) 

12 CONTINUE 

DO 13 1=1,6 

X TE Mp ( I )=X10( I ) -X 1 2 ( I ) 

13 CONTINUE 

CALL SUB1 (XTEMP,M2»Z4) 

DO 14 J = 1 , 6 

X 10 ( J + 6 ) =X10( J + 6) +XDEL( J) 

14 CONTINUE 

DO 15 J= 1 , 6 

XTEMPI J ) = X 1 0 ( J)-X13( J) 

15 CONTINUE 

CALL SUB1(XTEMP,M3,Z5) 

DO 16 J = 1 , 6 
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X10(J+6)=X10(J+6)+X0EL(J) 

16 CONTINUE 

00 1.7 J= 1 1 6 

XTEMPI J)=X13( J > -X 1 2 ( J) 

17 CONTINUE 

CALL SLB1(XTEMP,K23,Z6) 

DO ?9 J=l,6 

X12(J)=X12(J+6)-K3F23*XCEL{J) 

X13(J)=X13( J+6)+M2F23»XDEL( J) 

X10( J)=X10( J+6) 

29 CONTINUE 

DO 18 J = 1 » 3 
Y1C(J)=X1C< J)*CML 
Y12(J)=X12( J)*CML 
Y13(J)=X13( J)*CML 
Y10IJ+3) =X10( J+3) *CMV 
Y12< J+3)=X12< J+3)*CMV 
Y13 ( J + 3 ) =X 13 ( J + 3 > »CMV 

18 CONTINUE 

R10E2=X1C(1) »*2+X10 ( 2 ) *»2+XlC(3 ) «*2 

R20E2=(X1C( 1 ) -X12 ( 1 ) )«*2+(X10(2)-X12(2) )**2 + (X10(3)-X12(3) )**2 
R10=DSQRT ( R 10 E2 ) 

R20=DSGRT ( R20E2 ) 

Q=( M 2/tR10E2*R20E2) ) * ( M 1 /R 10+M 1 / R20 ) + 

1 (R2/36CC.C00)»(M1/R10*»3+M1/R2C**3)+ 

2 (M3/234CO.DOO** f 3)*(M2/R20E2+Ml/R10E2) 

R 1 0 = R 1 0*CM L 

R20=R20*CML 

1 F ( rpR.EC.OlGU TO b-j 

WRITE ( 3, 2C ) ( Y10I I ) , 1 = 1 .6) , TIMtD 
WRITE ( 3, 19 ) ( Y 1 2 ( I ) • 1 = 1 .6 ) 

WRI T E( 3, 19 M Y 1 3 ( I ) , 1 = 1.6) 

19 FORMAT ( IF ,6016.8) 

20 FOR M A T ( 1FC,7D16.8 ) 

WRITE (3, 28) Q, RIO, R20, TIRED 
28 FORMAT ( 1H ,4016.8) 

65 CONTINUE 

I F ( T I MED .LE.TIMEMX1GC TC 9 
STOP 

60 TIME D= TIMED N 
I PR =0 
GO TO 61 
END 

$ I B FTC SUFI LIST, REF .NOCECK 
SUBROUTINE SUB1 (X,M,Z) 

DOUBLE PRECISION 
D X ( 6 ) , M , Z , 

DX10,X12,X13,XTEMP,XKEP,XDEL, 

DM1, M 2,M3,M12,M13,M23,M2F12,M2F23,M3F13,M3F23, 

DCML,CMV,CMT,Y10,Y12,Y13,TIMED.TIMEMX, 

0TAUn,TAU,R10,R20,TEMP,Q,RlCE2,R2CE2, 

DZ1,Z2,Z3,Z4,Z5,Z6, 

DL,Cl,C2,C3,CELTA,H, 

DRE2,R,KS I ,ETA,ZETA,CHI 

COMMON X1C(12),X12(12),X13(12),XTEMP16) ,XKEP16),XDEL(6), 
CM1, m 2,M3,M12,M13,M23, M2F12, M2F23, M3F13,M3F23, 

CCML,CMV, CRT, Y10(6),Y12(6),Y13(6), TIMED, TIMEMX, 

CTAUC,TAU, RIO, R20, TEMP, C,R10E2,R2CE2, 

CZ1,72,Z3,Z4,Z5,Z6, 

CL»C1,C2»C3*DELTA»H» 
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CRE2,R,KSI,ETA,ZETA,CHI ,DUM< 1C00) 
RE2=X(1)**2+X(2)**2+X(3)**2 
R=DSQRT ( RE2 ) 

KSI=N*(TAU**2)/R**3 

ETA=(X(1)*X(A)+X(2)*X(5)+X(3)*X<6) )*TAU/RE2 
ZETA=IX(A)**2+X(5)**2+X(6)**2)*( TAU**2 ) /RE2-KS I 
CHI=KSI-ZETA 
CALL SUB2 (Z) 

DO AO 1=1,3 

XDEL II) = t-KSI»Z»*2) « (C2*X( I ) +C3* Z*TAU*X I 1+3 ) ) 

XDELI 1 + 3 ) = I -K S I *Z ) * I C 1*X ( I ) +C2 * Z «T AU* X ( 1 + 3) ) / (DELTA* TAU) 
XKEPt I )=X( I )+TAU*X( I+3J+XDEL ( I ) 

XKEPl 1+3 )=X( I +3 )+ XDEL l 1+3) 

AO CONTINUE 
RETURN 
END 

SIBFTC SUR2 L I ST , REF , NODECK 

SUBROUTINE SUB2IZ) 

DOUBLE PRECISION 
DZ, 

DX10,X12,X13,XTEMP,XKEP,XDEL, 

DM 1,^2, M3 ,f'12,M13,M23,M2F12,M2F23,M3F13,M3F23, 
DCML»CMV»CMT ,Y10,Y12,Y13,TIMEC, TI B EM X , 
DTAUC,TAU,R10,R20,TEMP,Q,R10E2,X2CE2, 

DZ1, 72,Z3,Z4,Z5,Z6, 

DL, Cl, C2,C3, DELTA, H, 

CRE2,R,KSI,ETA,ZETA,CHI 

COMMON X10(12),XI2(12) ,X13( 12) ,XTEMPI6) ,XKEP(6) , XDEL I 6 ) , 
CM1,M2,M3,M12,M13,M23, M2F12, M2F23, M3F13.M3F23, 

CCML ,CMV,CMT ,Y10( 6 ),Y12(6),Y13(b) .TIMED, TIME MX, 
CTAUD*TAU»R10*R20,TEKP,Q,R1C£2,R2CE2, 

CZ1,Z2,Z3,ZA,Z5,Z6, 

CL,Cl,C2,C3,0tLTA,H, 

CRE2,R,KSI ,ETA,ZETA,CHI ,DUM( 1CCO) 

I TER=0 

30 ITER=ITER+1 
L=CHI*Z**2 

C2=.5DOO* (M1-IL/12.D00) * ( M 1- I L/ 3 C . DOO ) * 

1 ( *1- (L/56.DOO) *( Ml- ( L/90. DOO ) ) ) ) ) 

C 3= l M1/6.CC0 ) *(Kl-( L/2C .000 ) * I N 1- l L/A2 .DOO ) * 

1 (M1-(L/72.DC0)*(M1-(L/11C.DCC) ) ) ) ) 

C1= M 1-L*C3 

DELTA=M1+C1*ETA*Z+C2*ZETA*Z**2 

H=C2*ETA*Z**2+C3*ZETA*Z**3+Z-M1 

Z = Z-H/CE LT A 

IFICABSIH) .L£. 1.0-0 7 ) RE TURN 
I F ( TTER.LE. 10)GC TC 30 
RETURN 
END 


$ I BS YS 

*JOB 05 1 C STUMPFF3 RH HILLIARC ATT. EHW RIGGS BLOG X-7267 

$CATE 102766 

SEXECUTE IBJOB 

$1 BJOB C-0, MAP, LOGIC, SCtRCE 

$ I BFTC MAIN LIST, REF, NOCECK 
COUPLE PRECISION 
DERR.TIMEDN, I TDN , I TO , A ,C X 10 , 

DX10.X12, X13,XTEMP,XKEP,XCEL, 
DST2X10,ST2X12,ST2X13,ST2X20,ST2X3C,ST2X23, 

CMi,M2,tf 3,M12,M13,M23,M2F12,M2F23 ,M3F13, M3F23, 

DCML , CMV, CRT, Y10,Y12,Y13, TIMED, TIMEMX, 

OTA UC,TAU, RIO, R20, TEMP, C,R10E2,R2CE2, 

DL,CI ,C2,C3,CELTA,H,XFR1C .XAFR10, 

OXFRI2,XFR13,XFR20,XFR23,XFR30,XAFR12,XAFRI3,XAFR20,XAFR23,XAFR30, 
0RE2,R,KSI ,ETA,ZETA,CHI , R i 2 , R 1 3 , R 23 , R30 , Z, 
0R01C,RC12,RD13,RCC10,RCC12,RC0I3, 

DR10F2C,R2CE2C,R10C,R20C 
DOU B L E PRECISION h Y, 10 , In V* 1 2 , R V* 1 3 , V R20 , W W 30 , WH 2 3 , 
DFF1C,FF12»FF13»FF20,FF3C,FF23 
COUPLE PRECISION C 
COUPLE PRECISION 

DW.TCTOAY, TDAYCU,ONED AY, SUT,TAUN, COUNT, AT,BT , 

CDEL10,CEL12,DEL13,RCNE,X20,DEL20,X30,DEL30,X23,DEL23, 

CEPS1C,XASTIC,EPSI2,XAST12,EPS13,XAST13*XAST2C,XAST23,XAST30, 

CRAST10,R<ST12,RAST13,RAST20,RAST23,RAST30,S10,S12,S13,S20, 

DS23,S3C,RCCT10,RDOT12,RCCT13,RC1CF,RC12F,RC13F,RCC10F,RCOI2F, 

0RC0l3F,RWGT,R0KGT,TSC240,R4FlO,R4F12,R4F13,TAU45,RC4FlO, 

CBX1C.BX12 ,BX13, 

DRD4F12,RC4F13,X4F10,X4F12,XAF13 
DOUPLE PRECISION E PS 2 0 , E P S 30 , EP S 2 3 

CCRRON X1C(12),X12(12),X13(12),XTEMP(6),XKEP(6),XCEL(6), 
CM1,M2,M3,R12,R13,M23, R2F12, M2F23, Rl 3F 1 3 , M 3F 2 3 , 

CST2XI0I6 ) , ST2XI2I6 ) ,ST2XI3 ( 6) , ST 2X201 6) ,ST2X3C< 6 ) ,ST2X23(6) , 
CCML,CRV,CRT,Y10(6),Y12(6),Y13(6),TIME0,TIMEMX, 

CTAUP ,TAU,R1C,R20,TEMP,C,R1CE2,R2CE2, 
CL,C1,C2,C3,CELTA,F,XFR1C(3),XAFR1C(3), 
CRE2,R,KSI,ETA,ZETA,CHI,R12(6),R13(6),R23(6),R30(6),Z, 
CR10F2C,R2CE2C,R10C,R20C 

COMMON V, Vm 1 0 t 6 ) » Vi I* 1 2 ( 6 ) » V* I* 1 3 ( 6) , In In 2 0 ( 6 ) ,kW30(6) ,WW23 ( 6 ) , 

CFFI0I6) ,FF12(6) ,FFI3(6) ,FF20(6) ,FF30(6) ,FF23(6) 

COMMON 

CVi,TOTCAY,TDAYCU,ONEDAY,SViT,TAUN, COUNT, AT, BT, 
CTSC’A0,BX10(6),BXI2(6),BX13(6), 

CRCC!2F{6)*RC013F(6) * R ImG T (A 1 ) , ROWG T ( A ) , R4 F 10 ( 3 ) ,R4F12f 3 ) , 

CRD 1 r ( 3 ) , RC 12(3), RD13 (3), RCD1C(3),RC 012(3), R0C13 ( 3) , 

CDEL10I6) , C E L 1 2 ( 6 ) ,DEL13 (6) ,CEL20 (6) ,0EL23( 6 ) , CEL 301 6 ) , 

CR0N B »X20(6 ) *X23(6) » X 3 0 ( 6 ) ,EPS10( 6),EPS12(6) ,EPS13(6), 
CXAST1C(6)»XAST12(6)*XAST13(6)»XAST20(6),XAST23(6),XAST30(6)» 
CRAST1C{6),RAST12(6),RAST13(6),RAST2C(6) , R AS T2 3 ( 6 ) , R AST30 ( 6 ) , 
CS10(6),S12(6),S13(6),S2C(6),S23(6),S30(6),R00T10(6), 
CRCCT12(6),RCGT13(6) ,RC1CF(6),R012F(6),RC13F(6),RCC1CF(6), 

CRAF 13(3) ,TAUA5*R0AF1C(3) , R0AF12I 3) , RDAF 13 ( 3 ) , XAF 10 ( 6 ) , XAF 12 ( 6 ) , 

C X AF 1 3 ( 6 ) ,XFR12(3) , X F R 1 3 ( 3 ) ,XFR20 (3),XFR23(3),XFR30(3),XAFR12(3), 
CXAFR13(3),XAFR20(3),XAFR23(3),XAFR3C(3) 

COMMON EFS20(6),EPS3C(6) »EPS23(6) 

C STUMPFF3 eY EH VstISS AND RH HILLIARC RIGGS BLCG 

900 REAP(2,8CC ) { Y 1 0 ( I ) , 1=1, 6 ) 

REAC(2,8CC) ( Y 1 2 ( I ) ,1 = 1,6) 

REAC(2»8CC ) ( Y 1 3 ( I ) ,1 = 1,6) 


22 



non 


800 FORMAT ( 3C16.8) 

REAP ( 2,8 C2)AT,BT,CNECAY,TIMEMX,MCRE 

802 FORMAT (4C16.8t7X, II) 

REAP ( 2 , 803 ) C 

803 FORMAT (DIE. 8) 

CML=6378. 16 5D00 

CMT = { ( 6 3 78, 165CCC**3 >/. 3986032DC6) 

CMT = D SCR T ( CMT ) 

CMIV = CML/CMT 

INPUT IS IN KILOMETERS AND SECONDS. 

MULTIPLY CANCMCAL BY C M L , CMT CR CMV, TC GBTAIN METRIC 
TAUT IS IN DAYS. TAU IS CANCNICAL. 

MONF= 1 . DCC 
M1=1.CC0 

M2= 1 .DCO/ei. 3C15CC52DC0 

M3=732951.2658DC0 

M12 = M 1 + M2 

M 1 3 = M 1 + M 3 

M23=M2+M3 

M2F12=M2/MI2 

M2F 23 = M 2 /M 2 3 

M3F I 3 = M 3 /M 1 3 

M3F?3=M3/M23 

TIMFD =o.cco 

COUNT=C. CCO 

TAU\=G.DCC 

0NECAY=(CNECAY*864CG.CCQ ) /CMT 
TOTE AY = 0.CC0 
S W T = C • DO C 

kRI Tt ( 3, 26 ) { Y10 ( I ) , 1 = 1,6 ) ,TIMED 
HRITE ( 3 t 19 ) ( YL2< I ) , 1 = 1 ,6) 

V*RITE< 3, 19 ) ( Y 1 3 { I ) , 1 = 1,6) 

26 FORMAT ( 1F1»6D16.8,D16.8) 

WRITEOt 8C1 )AT,Br,C, Ml 

801 FORMAT ( 4F A T = , D 1 6 . 8 , 2 X , 4 H B T = f D 1 6 . 8 , 2 X , 3H C=,C16.8,2X, 
15H m L = f c 1 6 • 8 ) 

RWG T ( 1 ) = 72 1 • DCO 
RWGT ( 2 ) =476. DOO 
RWGT ( 3 ) = 24 5 . DC 0 
R1«GT (4)=18.CC0 
RDhGT ( 1 ) =6 A . DCO 
RDtoGT ( 2 ) =24. DCO 
RDWM(3)=64.DC0 
RDtrfGT (4) =14. DCO 
W = 6 . DC C 
TAU=(^+AT)/BT 
9 CONTINUE 

C0UNT=C0UNT+1 .000 
TSQ 240= ( T A U * * 2 )/240.CC0 
TAU45=TAU/45.DC0 
DO 160 1=1,6 
RD 1 CF ( I )=C.CCO 
RD12F ( I ) =C • DC 0 
RD13F(I)=C.CC0 
RCDIOFCI ) = C . DOO 
ROD 1 2F ( I ) = 0 . DOO 
RCD 1 3F ( I ) = C . DCO 
350 CONTINUE 

CO 341 J = 1 » 4 
TAUN=TAUN+TAU 



DO ? 5 1 = 1,3 
BX 10(1) =Y 10(1 J/C^L 
BX12( I ) - V 1 2 ( I ) /CML 
0X13(1 ) = Y 1 3 ( I )/Ct*L 
8X 1 0 ( 1 + 3 ) = Y 10 ( I + 3 ) /CN V 
BX 1 2 ( I +3 ) = Y 1 2 ( r+3)/C^V 
BX13( I+3)=Y13( I+3)/CPV 
25 CONTINUE 
A = J 

TAU=A*TAU 
I PR = 1 

61 CONTINUE 
C 

C INPLT/m,X10 CUT PL T / X10,CEL10 

C 

CALL SLB1(BX10,H1) 

DO ICC 1=1,6 
X10U ) = X K E P ( I ) 

CEL10 ( I ) =XDEL ( I ) 

100 CONTINUE 
C 

C INPLT/yi+N2,X12 OUTPUT/ X12,DEL12 

C 

CALL SUB1 (BX12»K12) 

DO 11C 1=1,6 
X12(I)=XKEP(I ) 

DEL 12 ( I ) =XCEL ( I ) 

110 CONTINUE 
C 

C INPLT/P1+N3,X13 OUTPUT/ X13,CEL13 

C 

CALL SUBl(BX13,yi3) 

DO 120 1=1,6 
X 1 3 ( I ) = XKE P ( I ) 

D E L 1 3 ( I ) =XCEL ( I ) 

120 CONTINUE 

DO 13 1=1,6 

XTEVP( I ) =8X10 ( I ) - B X 1 2 ( I ) 

13 CONTINUE 

C 

C INPLT/P2, X1C-X1 2 OUTPUT/ X20,DEL20 

C 

CALL SUB1(XTEVP,^2) 

CO 130 1=1,6 
X 20 ( I ) = XK E P ( I ) 

DE L 20 ( I ) = X D E L ( I ) 

130 CONTINUE 

DO 15 P=l,6 

XTEVP(M) =BX10 (MJ-BX13 (fM 
15 CONTINUE 

C 

C INPUT/N3,X1C-X13 OUTPUT/ X30,DEL30 

C 

CALL SUB1 (XTEPP,M3) 

DO 1 4 C 1 = 1,6 
X30 ( I ) = X K E P ( I ) 

D EL 30 ( I ) =XC EL ( I ) 

140 CONTINUE 

DO 17 N = 1 , 6 

XTE*P(M)=BX13(M )-BX12(M) 


oooo oooo ooo ooo oooo ooo 


I 


17 CONTINUE 

I NPL T /M2 +M 3, X 1 3- X 1 2CUT PUT/ X23,DEL23 

CALL SLB1 < XTEMP.M23) 

00 150 1=1,6 
X23 ( I )=XKEP(I ) 

DEL 23 ( I ) = XD E L ( I ) 

150 CONTINUE 

CCMPUTE EPS10 E 1C = CEL2C+DEL 3C+M2/ ( M 1 + M2 ) «CEL12+ (M3/ ( M1+M3 ) 

•DELU3 

DO 190 1=1,6 

EPS 10(1 )=CEL2C( I ) +DEL 30 ( I > +M2F 1 2 «0EL 1 2 ( I) +M 3F 1 3*CEL 1 3 < I ) 

XAST101I ) = X 10 ( I ) +EPS 1C ( I ) 

COMPUTE EPSL2 E12=(M3/(M1+M3) >*DEL13-(M3/(M2+M3) )«DEL23 

EPS 12 ( I) =M3F13*DEL13( I ) -M 3F 2 3*0E L2 3 ( I ) 

XAS T 1 2 ( I ) = X 12 ( I ) + EPS 1 2 ( I ) 

CCMPLTE EPS03 t 1 3= ( M 2 / l M2+ M 3 ) ) *DEL 23+ I M2/ ( M 1 +M2 ) ) *DEL 1 2 

EPS13I I )=M2F23*DEL23( I ) + M2F 1 2*DE L 1 2 ( I ) 

XAST131 I ) = X 1 3 ( I ) + fc P S 1 3 ( I ) 

EPS20( I ) =E P S 1 0 ( I 1-EPS12I I ) 

EPS30( I )=EPS10( I )-EPS13( I ) 

EPS23( I )=EPS13( I )-EPS12( I > 

XAST20I I ) = XAST1C( I> — XAST121 I ) 

XAST3C ( I >=XAST10( I J-XAST 1 3 ( I ) 

XAST2 3II ) = XAST13( I )-XAST12( I ) 

190 CONTINUE 

COMPUTE R ( I J ) SCUAREC ANC RAST(IJ) SOUAREO 

R I I J ) «*2= X< I JM 1>**2 + X( I J ) (2)**2 + X( IJ) ( 3>**2 

R10 = DSCRT(X10(1)*X10(1) + X1C(2)*X1C(2)+X10(3)-»X1C(3)) 
R12=DSCRT(X12(1)*X12(1)+X12(2)*X12(2)+X12(3)»X12(3)) 
R13=DSCKT(X13(1)*X13(1)+X13(2)*X13(2)+X13(3)*X13(3)1 
R20=DSCRT(X20(1)*X20(1)+X2G(2)»X2C(2)+X2G(3)*X20(3)) 
R23=CSCRT(X23(l)*X23(l)+X23(2)*X23(2)+X23<3)*X23(3)) 
R30=CSCRT(X3C(1)*X30(1)+X3C(2)*X3C(2)+X30(3)*X3C(3)) 

RAST10 = DSCRT ( XAST10I 1)«*2+XAST10 (2) «*2 + XAST10( 3 ) * »2 I 
RAST12=DSCRT(XAST12(1)*»2+XAST12(2)«*2+XAST12<3)*»2) 
RAST13=DSCRT(XAST13(1)»*2+XAST13(2)«*2+XAST13(3)«»2) 
RAST20=DSCRT(XAST2C(1)**2+XAST20(2)«*2+XAST20(3)«*2) 

RAST23 = OSCRT ( XAST23 ( 1 ) ««2 + XAST23 (2 ) ««2+XAST 23 ( 3) « » 2 ) 
RAST3C=DSCRT(XAST30(1)«*2+XAST30(2)«»2+XAST30(3)»»2) 

COMPUTE S10,S12,S13,S20,S23,S30 V.FERE IN GENERAL 

S(IJ)=X(IJ)/R(IJ)»*3*XAST(IJ)/RAST(IJ)**3 


DO 250 1=1,3 

SIO ( I )=X1C ( I ) /(R1C*R1C*R10)-XAST1C( I )/RAST10*«3 
S12(I)=X12(I)/(R12*R12*R12)-XAST12(I)/RAST12*«3 
S13(I)=X13(I)/(R13«R13»R13)-XAST13(I)/RAST13«*3 
S20(I)=X2C(I)/(R20*R2C»R20)-XAST20(I)/RAST20*«3 
S 2 3 ( I ) = X23 ( I ) /(R23*R23*R23)-XAST23( I )/RAST23*«3 
S 30 m = X3Cm/(R3C*R3C*R2C)-XAST3C(I)/RAST30**3 


25 


I 


oooooooon 


CONTINLE 


250 


COMPUTE ROOT (IJ) 

ROOT I I J » = C M i I) + M( J) ) «S(I J) ♦f'(K)«(S( IK)+S(KJ> ) 
+K(L )*(S( IL)+S(LJ) ) 
lnHERE I .NE. J .NE. K .NE. L 
AND (I,J,K,L)=(C,1,2,3) 


DO 260 1=1*3 

ROOT 10 { I ) = K1*S10( I)+K2»(S12( I ) + S20( I ) )+K3*(S13 (I) + S30(I) ) 

ROOT 1 2 ( I ) = M2*S12( I )+R3» (S13( I ) — S 23 t I ) ) 

RDCT13I I ) = f*13*S13( I > + R2*(S12(I) + S23( 1)) 

260 CONTINLE 
TAL=TAL/A 
DO 340 1=1,3 

RD1C ( l ) = RViGT(J)*RCCTlC( I ) 

RD1CFI I )=RD10F( I ) + RO 1 0 ( I ) 

RD1 2 ( 1 ) = hV<GT(J)*RD0T12( l ) 

RD12F1 I J = P C 12 F ( I ) + R U 1 2 ( I ) 

RD13( I ) = RViGT(J)*RDCT13( I ) 

RO 1 3F ( I )=RD13F( I 1+RD13I l ) 

RDD 1 0 ( I ) =RDWGT ( J ) «RDCT1C ( I ) 

RCDlOFtl )=RCD10F( I )+RCDlC( I ) 

RCD12 ( I ) = R D in G T ( J ) *RDOT 1 2 ( I ) 

RCO 1 2 F ( I )=RCC12F( I )+RCD12( I ) 

RDD 1 3 ( I )=RDHGT ( J ) *RCCT13 ( I ) 

RCD13F ( I )=RCD13F ( I J+RCD 1 3 ( I ) 

340 CONTINUE 

341 CONTINUE 

DC 370 1=1,3 

R4F10 ( I )=TSQ240*RC1QF( I ) 

R4F12I I ) =TSG240*R012F ( I ) 

R4F131 I )=TSG240*RC13F( I ) 

RD4F1CI I ) = TAU45*RCD1CF( I ) 

RD4F12I I ) = TAU45*RCD12F t l ) 

RD4F13I I ) = TAU4 5*RCD13F ( I ) 

370 CONTINUE 

CX IOC SORT ( (R4F 1C ( 1 ) * «2 ) + ( R4F 1C { 2 ) * *2 ) + ( R4F 10 ( 3) **2 ) ) 

380 00 390 1=1,3 

X4F 1 0 ( I ) =XA S T 10 ( I ) +R4F1C ( I ) 

X4F12I I ) =XAST12( I ) + R4 F 12 ( I ) 

X4F13II ) = XAST13( I ) +R4F 1 3 ( I ) 

X4F101 I+3)=XAST10(I+3)+RC4F10(I! 

X4F12I I + 3)=XAST12 ( I + 3)+RC4F12( I ) 

X4F13I I+3)=XAST13 ( I + 3 H-RC4F 1 3 < I > 

390 CONTINUE 
69 DG 18 J= 1 , 3 

Y1C(J)=X4F10(J)*CNL 
Y12(J)=X4F12( J)*CRL 
Y 1 3 f J) = X4F13( J)*CNL 
Y10(J+3)=X4F1C(J+3)*CNV 
Y12 (J+3)=X4F12( J+3)*CRV 
Y13 ( J+3) =X4F13( J+3) *C N V 
18 CONTINUE 

R10F2=Y1C(1)*»2+Y1C(2)*»2+Y1C(3)«»2 

R2CF2=(Y1C(1)-Y12(1))«*2+(Y1C(2)-Y12(2) )**2+(Y10(3)-Y12(3))**2 
R10=CSCRT(R10E2) 

R20=DSCRT (R20E2) 



R1CE2C=X4F10( 1)»»2+X4F10(2)*«2+X4F10(3) « * 2 

R20P2C = ( X4F 10 ( 1 1-X4F12 ( l))«*2+(X4F10!2)-X4F12!2) )«*2+ (X4F 10 ( 31-XAF 

112 m )»*2 

R10C=DSQRT (R1CE2C ) 

R20C=DSQRT (R20E2C ) 

W=C w lNl(R10C,C»R20C) 

TAU'=!W + AT)/BT 

IFt ShT-1 .CCC) 1CC6.1C01. icoe 
1006 I F ( PNEOAY- ( T AUN + A » T AL ) )1CCC, 1005,9 
1005 SWT=i.COC 
GC TO 9 

10C0 TAU= !ONEOAY-TAUN ) /A 
SWT = 1 .000 
GO TO 9 

1001 TOTCAY=TCTDAY+ONECAY 
TDAYCU=(TCTCAY»CMT)/864CC.CCC 
SWT = 0 .COC 

T I N EC = TO A YCU 

IF( I PR. EG. 0) GC TO 65 

WRITE ( 3,2C M Y10( ! ) , 1 = 1,6) ,TINED 

WRITE(3, 19) { Y 1 2 ( I) ,1=1,6) 

WRITE! 3, 15 ) (Y13< I) , 1=1,6) 

19 FORMAT! IF ,6016.8) 

20 FGRFAT < 1FC , 7016.8 ) 

2050 FORMAT ( IF 3C2A.16) 

WRITE (3, 26 )CX10,R1C,R20,TIFEC,C0UNT 
28 FORMAT ( IF ,5016.8) 

IF! ICINT (TIPEKX)-IDINT ( TDAYCL) ) 1003,1003,1002 

1002 CONTINUE 

T AUN = 0 . DC C 
GO TO 9 

1003 IF! WORE. EC. 1 )G0 TC 9C0 
65 CONTINUE 

STOP 

END 

SIBFTC SUP1 LIST, REF, NCCECK 

SUBROUTINE SUB1 ! X , N ) 

COUPLE PRECISION 
OX ( 6 ) « P » 

CERR.TIPECN, ITDN, ITC,A,CX10, 

DX1C,X12,X13,XTEPP,XKEP,XCEL, 

CST2X10,ST2X12,ST2X13,ST2X20,ST2X2C,STZX23, 

DP1,P2,M3,P12,P13,P23,P2F12,F2F23,P3F13,P3F23, 

0CML*CMV»CPT*Y10,Y12»Y13,TINEC»TIPEPX, 

CTALC,TAU,R1C,R20,TEMP,C*R1CE2,R2CE2, 

DL.Cl ,C2,C3,CELTA,H,XFR1C,XAFR10, 

DXFR12,XFR13,XFR20,XFR23,XFR3C,XAFR12,XAFR13,XAFR20,XAFR23,XAFR30, 

DKE2,R,KSI,ETA,ZETA,CFI,R12,R13,R23,R30,2, 

DR10F2C,R20£2C,R10C,R2CC 
COUPLE PRECISION W W 1C , W W 1 2 , W W 13 , V W2C , WW30 , W W2 3 , 
DFF10,FF12,FF13,FF20,FF3C,FF23 
OOUPLE PRECISION 

DW,TCTOAY,TCAYCU,ONEDAY,SWT,TAUN,CCUNT,AT,BT, 

DRD1C,RD12,R013,RCC10,RDC12,RCD13, 

CDELlO,0ELl2*DEL13*F l ONE*X2O»CEL2C»X3C»DEL3O»X23»DEL23» 

DEPS10,XAST10,EPS12,XAST12,EPS13,XAST13,XAST20,XAST23,XAST30, 

DRAST10,RAST12,RAST13,RAST2C,RAST23,RAST30,S10,S12,S13,S20, 

DS23 ,S30, RCOT10,RDOT12,RDCT13,RD1CF,RD12F,RD13F , RCD1 OF , RDD12F , 
DRCD13F,RWGT , RCWGT ,TSC24C,R4F10,RAF12,R4F13.TAU45,R0AF10, 
DBX1C,BX12,BX13, 



0RD4F12,RD4F13,X4F10,X4F12,X4F13 
DOUeiG PRECISION EPS2C,EPS30,EPS23 

CORDON X1C(12),X12(12)»X13(12)»XTERP(6) ,XKEP 1 6 ) » XDEL ( 6 ) , 
CR1,R2,M3,R12,R13,R23, R2F12, R2F23, M3F 13 , R 3F 2 3 , 
CST2X10(6),ST2X12(6),ST2X13<6>,ST2X20<6),ST2X3C(6),ST2X23(6), 
CCML,CRV,CRT,Y10(6) ,Y12(6),YI3(6) »TIRED» TIREMX , 

CTAUC,TAU,R 10.R20, TERP , C , R 10E2 , R2 CE2 , 
CL,C1,C2,C3,DELTA,H,XFR1C(3),XAFR 10(3) , 

CRE2,R,KSI ,ETA,ZETA,CHI f R12I6) ,RI3(6) ,R23< 6) ,R30< 6) ,Z f 
CR10E2C,R2CE2C ,R10C,R20C 

COR PUN WW1C (6) ,V»W12(6) ,hhl3(6) ,hk2C(6) ,V>V»3C(6) ,Wh23<6) , 

CFFIC(6) *FFI2(6) »FFI3(6) , FF2C ( 6 ) , FF30 ( 6 ) ,FF23(6> 

CCRRON 

CVi»TCTCAY»TCAYCU*ONEOAY,SWT ,TAUN»CCUNT |AT »BT , 

CTSC 240,6X10(6) ,BX12( 6 ) ,6X13(6) , 

CRCD12FI 6) ,RCD13F( 6) ,RV<GT(4) ,RDWGT(4) ,R4F10( 3) ,R4F12( 3) , 

CRD10(3) ,RC12(3) ,RC13(3) ,RCC10(3) , RCD 12 ( 3 ) »RCC13(3), 

COEL 1 0(6) , C E L 1 2 ( 6) , DEL 1 3 (6) , DE L20 ( 6 ) , DEL 2 3 ( 6 ) , CEL 3C ( 6 ) , 

CRONF , X20 ( 6 ) ,X23(6),X30(6),EPS10(6),tPS12(6),EPS13(6), 
CXAST10(6),XAST12(6),XAST13(6),XAST20(6),XAST23(6),XAST30(6), 
CRAST10(6),RAST12(6),RAST13(6),RAST20(6) , R A S T2 3 ( 6 ) , R A S T30 ( 6 ) , 

CS10 (6) , S 12 ( 6 ) , S 1 3 ( 6 ) , S2 0 ( 6 ) , S 23 ( 6) , S3C ( 6) ,RDCT10(6) , 

CRD0T12 (6 ),ROOT13( 6) , RD 1CF ( 6 ) , R01 2F ( 6 ) , R C 1 3F ( 6 ) , RCC 10F ( 6 ) , 

CR4F 13(3) ,TAU45,R04F10( 3 ) , RC4F 12 ( 3 ) , RD4F 13 ( 3 ) ,X4F1C(6) ,X4F12 I 6) , 

C X4F 1 3 ( 6 ) , XFR12 ( 3) »XFR13( 3) ,XFR20 (3) ,XFR23(3) ,XFR3C ( 3) ,XAFR12(3) , 
CXAFR13(3) ,XAFR2C(3) ,XAFR23(3) ,XAFR30(3) 

CCR V QN EPS2C(6) ,EPS30(6) ,EPS23( 6 ) 

RE2 = X(l)«*2 + X(2)-**2 + X(3)**2 
R = OSGRT( RE2 ) 

KSI=M*(TAU«*2)/R»*3 

ETA=(X(l)«X(4)+X(?)*X(5)+X(3)*X(fc))»TAU/RE2 
ZETA=<X(4)«*2+X(5)»*2+X(6)«*2)»( TAU*»2 ) /RE2-KS I 
CHI=KSI-ZETA 

Z=MCNE-. 5CC0*ETA-{ l.OCC/ 6 . DOO ) * ZETA+ . 5 COO* ETA* *2+ ( 5 . 000/ 12 . DOO ) 
1 *ET A » Z E T A 
I TE R =0 

30 I TER=I TER + 1 
L =CH I » Z * *2 

C2=.5D00*(R0NE-(L/12.CCQ)*( RONE- (L/30.0C0)* 

1 < RONE- (L/56.D00)* (RONE- (L/90.CC0) ) )) ) 

C 3= (NONE/ 6.0C0) * (RONE-( L/ 20. DOO )« ( RONE- ( L / 42 . COO ) * 

1 (R0NE-(L/72.C00)*(RCNE-(L/110.CC0) ) ) ) ) 

C 1 = v ONE-L »C 3 

DELTA=R0NE+C1*ETA*Z+C2*ZETA*Z**2 

H=C?«ETA«Z«*2+C3*ZETA»Z»+3+Z-MON£ 

Z=Z-H/CELTA 

I F ( CABS ( H ) . LE . 1 . 0-07 ) GC TC 31 
I F ( ITER . LE . 10 ) GO TO 30 
31 CONTINUE 

DO 40 1=1,3 

XDEL ( I )=(-KSI*Z»*2)*(C2«X(I ) + C3« Z*TAU*X ( 1 + 3 ) ) 

XDEL( I+3)=(-KSI»Z ) *(C1*X ( I ) +C2*Z *T AU* X ( 1+3) ) / <OELTA*TAU) 

XKEP( I ) = X(I) + TAU*X( I + 3) + XDEL( l ) 

XKE P ( I+3)=X( I+3)+XDEL( 1+3 ) 

40 CONTINUE 
RETURN 
END 

$DATA 

.23621551006 -. 264078 1 5CC6 1 5944798CC6 

. 8441 5966DCC -. 71 504 167C-1 4078 1976C- 1 
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I 


• 2284 1 895006 

-.2861 2 319DC6 

-• 1 6078232C06 


•8243782 iOCC 

• 50 846563CC0 

. 18637202CC0 


-•31957216008 

• 1 364 2 1 3 1CC9 

.59 15 7903 COB 


-.28628154002 

-.56407282001 

-.2*463687001 


2.0DC0 

19.0CC0 

LO.OCCO 

179.0000 

l.CDCC 
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$ I B SY S 

$ JOB 0510 STUMPFF4 E.H. WEISS RIGGS BLOG X-7266 

SDATE 102766 

4EXECUTE I 8 JOB 

SIB JOB GO, MAP, LOGIC, SOURCE 

$ I BFTC MAIN LIST, REF, NGDECK 

DOUBLE PRECISION 
DERR.TIMECN, I T ON , I TO , A , CX 10 , 

DX10,X12, X13,XTEMP,XKEP,X0EL, 

DST2X10,ST2X12,ST2X13»ST2X20,ST2X3C»ST2X23, 
DM1,M2,M3,M12,M13,M23,M2F12,M2F23 , M3F 13 , M3F2 3 , 

DCML, CM V,CMT,Y10,Y12, Y 1 3 , T I MED , T I MEMX , 

DTAU0,TAU,R10,R20, TEMP, Q, RICE 2 ,R2CE2 , 

DL, C1,C2,C3,OELTA,H,XFR10,XAFR10, 

0XFR12,XFR13,XFR20 , XF R23 , XF R30 , X AFR 1 2 , X AFR 1 3 , X AFR20 , X AFR 2 3 , X AF R30, 
0RE2,R,KSI ,ETA,ZETA,CHI , R 12 , R 1 3 , R 23 , R30 , Z , 
ORD10,RD12,R013,RCCIO,RDD12,RC013, 

DR10F2C,R20E2C,R10C,R20C 
DOUBLE PRECISION W W 10 , WW 1 2 , W W 13 , W W2C , WW 30 , W W2 3 , 
DFFlf',FF12,FF13,FF20,FF3C,FF23 
DOUBLE PRECISION C 
DOUBLE PRECISION 

DW, TPTOAY, TDAYCU.ONEDAY, SWT , TAUN , COUNT , AT , BT , 
DDEL10,DEL12,D£L13,M(jNE,X 20,DEL20 ,X 30 , DE L30 , X2 3 , DEL2 3 , 
DEPS!0,XAST10,EPS12,XAST12,EPS13,XAST13,XAST20,XAST23,XAST30, 
DRAST10,RAST12,RAST13,RAST20,RAST23,RAST30,S10,S12,S13,S20, 
DS23,S30,RC0T10,RUCT12,RD0T13,RDICF,R012F,RD13F,RCC10F,RDD12F, 
DRDDl3F,RwGT,RDWGT ,TSC2 AO, RAF 10, RAF 12 , K4F 13 , T AU 45 , RD4F 10 , 

DBX10, 3X12,6X13, 

DDEL T , CEL (4 ) , 

DRD4F12,RD4F13,X4F10,X4F12,X4F13 
DOUBLE PRECISION T AUP , R 10F 4T , R 1 0 F T4 

COMMON X 1 0 I 12) ,X12( 12) , X 1 3 t 12! ,XTEMP(6) ,XKEP(6! ,XDEL(6) , 

CM1, V2.M3 ,M12,M13,M23, M2F12, N2F23, M3F13.M3F23, 
CST2Y10(6),ST2X12(6) ,ST2X13(6) ,ST?X20(6) .ST2X3016) ,ST2X23(6) , 
CCML,CMV,CMT , Y 10 I 6 ) , Y 12 ( 6 ) , Y 13 t 6 ) , TIMED, TIME MX, 

CTAUO, TAU.R10, R20, TEMP,C,R1CE 2..R2CE2, 

CL»C1,C2,C3, DELTA, H,XFR10(3) , X AFR 1C ( 3 ) , 

CRE2,R,KSI,ETA,ZETA,CHI,R12(6),K13(6),R23(6),R3C(6),2, 

CR10F2C,R2C£2C,R10C,K2GC 

COMMON WW10(6) , W W 1 2 ( 6 ) ,WW13(6) , W W20 < 6 ) * W w 30 ( 6 ) ,WW23(6) , 

CFF1C16) ,FF12(6),FF13(6) ,FF20(6),FF30(6) ,FF23(6) 

COMMON 

CW, TCTDAY.TCAYCU.ONEOAY, SWT , TAUN , CCUN T , AT , BT , 

CTSC?40,BX10(6) , BX12( 6) ,8X13(6) , 

CRCD12F(6),RCD13F(6),RwGT(4),RD«GT(4),R4FlC(3),R4F12(3), 

CRD 10(3), RC12( 3) ,RD13( 3) ,RCD10(3) ,RCD12I3), ROD 13(3), 

CDEL10(6) , CEL 1 2 ( 6 ) ,DEL13(6) ,DEL20 (6) , DEL 23 ( 6 ) , CEL 30 ( 6 ) , 

CMCNB , X20 (6 ) ,X23(6) ,X3C (6) ,EPS10( 6) ,EPS12(6) , E P S 1 3 ( 6 ) , 

CXAST10(6 ),XAST12!6! ,XAST13(6) ,XAST2C(o) ,XAST23(6) ,XAST30(6) , 
CRAST10(6 ) ,KAST12( 6) ,RAST13 (6) ,RAST20(6) ,RAST23 ( 6) ,RAST30(6) , 
CS10(6) , S 12 ( 6 ) , S 1 3 1 6 ) >S20(6) , S23 ( 6 ) ,S30(6) .RD0T1016) , 

CRD0T12 (6) ,RD0T13( 6) , RD1 CF ( 6 ) , RD1 2F ( 6 ) ,RD13F(6) ,RCD10F(6) , 

CK4F13 ( 3) ,TAU45,RD4F10( 3) ,RC4F12 (3) ,RD4F13( 3 ) , X4F 10 ( 6 ) , X4F 12 ( 6 ) , 
CX4F 13(6 ) , X F R 1 2 ( 3 ) , X F R 1 3 ( 3 ) ,XFR20(3) »XFR23(3) , X F R 3 0 ( 3 ) ,XAFR12( 3) , 
CXAFR13(3),XAFR20( 3 ) , XAFR23 ( 3 ) , X A FR30 ! 3 ) 

C STU V PFF 0 Y E.H. WEISS RIGGS BLDG 

900 REAr(2,8CC) (Y101I ) , I = 1 , 6 ) 

READ(2»8CC) ( V 1 2 ( I ) , 1 = 1 ,6 ) 

READ(2»8CC) (Y13II),I=1,6) 
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800 FORMAT ( 3C16.8) 

REAC( 2,802) A T , 8 T , CNt U A Y , T I ME M X , MCRE 

802 FOR M AT (4C16.6,7X, II) 

803 FORMA T ( D 1 6 . 8 ) 

CML=6378 .165000 

GMT = { (63 78. I 6 5000* *3)/. 3986 032006) 

CMT=DSGRT (CM ) 

CM V = CML / CM T 

I N D U T IS IN KILOMETERS AND SfcCCNDS. 

MULTIPLY CANONICAL bY C M L » CMT CR CMV, TC OBTAIN METRIC 
TAUC IS IN DAYS. TAU IS CANONICAL. 

I PR = 0 

MONE= 1 * DCO 
M1M.D00 

M2* l .DC0/S1.3C15OC52OC0 

M3 = 33295 I.26580C0 

M 1 2 = M 1 + M 2 

M I 3=M 1 +M 3 

M23=M2+M3 

M2F12=M2/M12 

M2F ?3 = M2 /M23 

M3F13=M3/MI3 

M3F?3=M3/M23 

TIMFD =O.CCC 

COUNT=O.CCO 

TAUN=C.OCC 

ONEDA Y = ( 0NEDAY*864CC.C00)/CMT 

TQTCAY=O.DCO 

SWT=0 . DOC 

WRITE (3, 26 ) { YlOU) , 1 = 1 ,6) , TIMED 
WRITE (3, 19) { Y 12 ( I ) , 1 = 1, 6) 

WRITE(3»19) ( Y 1 3 ( I ) f 1 = 1,6) 

26 FORMAT ( 1H1,6D16.8,D16.8) 

WRITEI3, 801 )AT,BT,C,M1 

80 1 FORMAT { 4HCAT=tD16.3,2X,4H B T = , 0 I 6 . 8 , 2 X , 3H C=tCl6.8 f 2X, 
15H Ml =,D16.8) 

RWGT( 1 ) =721.000 
RWGT ( 2)=476.D00 
RWGT { 3) =245.000 
RWGT ( 4 ) = 18 . COO 
ROWGT ( 1)=64.D00 
RDWGT ( 2 ) =2 A * DOO 
RDWGT ( 3) =6 A . 000 
ROWGT (4 ) =14.000 
DEL ( 1 )=-4.DC0 
DEL ( 2 ) = 6 • DC 0 
DEL ( 3 ) = — 4 . DC 0 
DEL (4 ) =1 .CCG 
W = 6 * 0 C 0 
TAU=BT 

9 CONTINUE 

C0UNT=C0LNT+1 .DCO 

TSC?40=< TAb**2> /240.DC0 

TAU45=TAU/45.D00 

00 '350 1 =1,6 

RD 1 OF ( I ) =C . DO 0 

RD 1 2F ( I )=C.CCO 

RD1 3F ( I )=C.COO 

RCD10F ( I ) =C • DCO 

RDD 12F t I ) = 0 . DGO 



ROD 1 3 F ( I ) = 0 . DOO 
350 CONTINUE 

DEL T = Q . DCC 
DO 34 L J = 1 * 4 
TAUN=TAUN+TAU 
DO 25 1=1,3 
BX1CC I ) = Y 1 0 ( I )/CVL 
B X 1 2 { I ) = Y 1 2 ( I )/CML 
BXI3( I ) = Y 1 3 ( I ) /CM L 
BX1C( I+3)=Y10(I+3)/CMV 
8X1?( I + 3 ) = Y 1 2 ( I + 3 ) /C M V 
B X 1 3 { 1 + 3 ) = Y 1 3 ( I + 3 ) /U M V 
25 CONTINUE 
A = J 

TAU = A*T At 
61 CONTINUE 
C 

C I NPUT/M 1 , X 1 0 ULTPUT/ X1G,DEL10 

C 

CALL SUB1(BX10,M1) 

DO 100 1=1,6 
X10 ( I ) = XKEP { I ) 

DEL10C I ) =XDEL ( I ) 

100 CONTINUE 
C 

C INPLT/M1+M2,X12 UUTPUT/ X12,DEL12 

C 

CALL SUB i(BX12»M12) 

DO 110 1=1,6 
X 12 ( I ) =X K E P ( I ) 

DEL12C I ) = XC EL ( I ) 

110 CONTINUE 
C 

C INPUT /Ml+M3fX13 OUTPUT/ X13,DEL13 

C 

CALL SUBl(6X13tM13) 

DO 120 1=1,6 
X 1 3 ( I ) =XKE P ( I ) 

DEL 1 3 ( I ) =XDEL { I ) 

120 CONTINUE 

DO 13 1=1,6 

XTE W P ( I ) =BX 1 0 ( I ) -BX 1 2 ( I ) 

13 CONTINUE 

C 

C INPUT/M2,X10-X12 OUTPUT/ X20,DEL20 

C 

CALL SU61(XTEMP,M2) 

DO 130 1=1,6 
X20U )=XKEP( I ) 

DEL 20 ( I ) =XD EL ( I ) 

130 CONTINUE 

DO 15 M = 1 , 6 

XTEMP(M) =exiO(M)-BX13( M) 

15 CONTINUE 

C 

C INPLT/M3,X1C-X13 CUTPUT/ X30,CEL30 

C 

CALL SU B 1 ( X TE MP , M 3 ) 

DO 140 1=1,6 
X30 { I ) =XKE P ( I ) 

DEL^O ( I ) =XCEL ( I ) 
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140 CONTINUE 

DO 17 M= 1 , 6 

XTEMP(M)=BX13(M)-BX12<M) 

17 CONTINUE 

INPUT/M2+M3.X13-X12ULTPLT/ X23.0EL23 

CALL SUB1(XTEMP,K23) 

DO 150 1=1,6 
X23 ( I ) =XKEP ( I ) 

DEL ?3 ( I ) =XCEL ( I ) 

150 CONTINUE 

C 0 M ° U T E EPS10 E10=CbL20+DEL30+M2/l MI+M2) *DE L 1 2+ l M3/ l M 1 +M3 1 

*DELU3 

DO 190 1=1,6 

EPS 10 (I ) =CEL20( I 1 +DEL30II) +M2F 1 2 «CEL 1 2 ( I ) +M3F 13*0EL 1 3 (I) 

XASTlOt I ) = X 1 0 t I ) + E PS 1 0 ( I ) 

COMMUTE FPSU2 E 1 2= ( M3 / ( M 1 + M 3 ) ) ‘DEL 1 3- ( M3/ ( M2 + M3 ) ) *DEL2 3 

E PS 1 2 ( I )=M3F13*DEL13( I )-M3F23*DEL23( I ) 

XAST121 I ) = X L 2 ( I )+EPS12( I ) 

COMPUTE EPSU3 t 1 3= < M2 /( M2+ M 3 >> *D£L 23+ ( M2/ I M1+M2 )) *DEL 12 

EPS131 I )=M2F23*DEL23( I )+M2F12*DELI2( I ) 

X A S T 1 3 ( I ) = X 13 (IJ + FPS 1 3 (I) 

XAST20II )=XAST10( I )— XAST12( I ) 

X A S T 3 0 ( I ) = XAST1C( I )— XAST 13( I ) 

X AS T2 3 ( I ) = XAST13( I ) - X A S T 1 2 { I ) 

190 CONTINUE 

COMPUTE R ( I J ) SGUARLC AND RAST(IJ) SQUARED 

R ( I J ) **2 = X( IJ) ( 1 ) ** 2 + X I IJ) (2)**2+X(IJ) (3)**2 

R10 = DSCRT< X10( 1)*X10 <11 +X10 ( 2 )*X I0( 2 ) +X10I 3) *X10 ( 3) ) 

R12=USQRI ( X12 (1)*X12(1)+X12(2)*X12(2)+X12(3)*X12(31) 

R13=DSCRT< X13I1)*X13 ( 1 1+X13 < 2 )*X 13 ( 2) +X13 t 3) *X13 ( 3) ) 
K20=DSQRT(X20<l>*X2u<l>+X2C(2)*X20<2)+X20(3)*X20(3)) 
R23=0SCRT<X23(l)*X23(l)+X23(2)*X23(2)+X2i(3)*X23(3)> 

R30 = 0SCRT ( X30 I 1 ) *X30( 1 ) + X3C ( 2 ) *X3C( 2 1 + X30 t 3 ) *X3C ( 3) ) 
RAST1C=DSCRT(XAST1C(1)**2+XAST10(2)**2+XAST10(3)**2) 
KAST12=DSCRT(XAST12(L)**2+XAST12I2)**2 +xAST12(3)**2> 

RAST13=DSCRT ( XAST 13 ( i )«*2+XAST13 (2)**2+XAST13(3)**2) 
RAST20=DSCRT(XAST20(1)**2+XAST20(2)**2+XAST20(3>**2) 

RAST23 = DSGRT ( XAST231 i ) **2+XAST23 (2) **2+XAST23( 31**2) 

RAST30 = 0SCRT< XAST 30 ( 1 ) **2+XAST30 (2 ) **2+XAST30 ( 3) «*2 ) 

COMPUTE S10,S12,S13,S20,S23,S30 *FIERE IN GENERAL 

S( I J 1 = X ( IJ)/R(IJ)**3*XAST( IJ) /RAST ( I J 1**3 


DO 250 1=1,3 

S10(I)=X1C(I )/(R10»k10»R10)-XAST1C1 I )/RAST10*«3 
S12 ( I 1 =X 12 ( I )/(R12*Kl2*R12)-XASTl2( I ) /RAST 12** 3 
S13( I 1 = X 1 3 ( I 1 /(R13*R13*R13)-XAST13( I )/RAST13**3 
S20II )=X2C( I 1 /(R20*R2C*R20 1-XAST2CI 1 1 /RAST20**3 
S23(I ) = X 2 3 ( I )/(R23*K23*R23)-XAST23( I)/RAST23**3 
S 30 m=X3C(I)/(R3C*l<30*R30)-XAST2C(I)/RAST 30**3 
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250 CONTINUE 


COMPUTE RCOT (IJ) 

ROOT ( I J) = (M( I )+M( J) ) *S( I J)+M(K)*< S( IK)+S(KJ) ) 
+M(L)*(S( ILJ+SILJ) ) 

INHERE I #NE . J .NE. K .NE. L 
AND (I,J,K,L>={0,1,2,3> 


DO 260 1=1,3 

ROOT 1 0 ( I ) =M 1* S 1 C ( I )+M2*(S12( I ) + S2Ct I ) >+N3*< SI 3 ( I ) + S30 ( I ) ) 

ROOT 1 2 ( I )=M12*S12 ( I ) +M3* ( S 1 3 ( I ) — S 2 3 ( I ) ) 

ROOT 1 3 ( I >=M13*S13(1)+N2*<S12( I )+S2 3(I) ) 

260 CONTINUE 

DEL T = RDOT 10 ( 1 ) **2 + RDOT 1 0 ( 2 )* *2+K CC T 10 t 3 )* *2 
DELT = DSORTICELT)*DEL(J)-*-CELT 
TAU=TAU/A 
DO 340 1=1,3 

RDlC(I)=RkGTlJ)*RD0T10<I) 

RD1CFI I ) =RC10F< I ) +RD10I I) 

RD12(I>=RNGT( J)*RDUri2(I) 

RD12FI I )=RD12F( I ) + RD 1 2 ( I ) 

RDl 3 ( I ) = RViGT(J)*RD0T13( I ) 

RD 1 3F ( I )=RC13F( I )+RC13( I ) 

RCD101 I ) =RDwGTI J ) *R0GT10 ( I ) 

RCD10FI I )=RCD10F( I )+KCDlC( I ) 

RDD12 ( I ) =RCWGT ( J ) *RDUT12( I ) 

RCD12F ( I ) =RCD12F ( I )+RCC 12J I ) 

ROD 1 3 ( I ) =R DWG T ( J ) « ROOT 13 ( I ) 

RDD13FII ) = RC013Ft D+RCD13 1I ) 

340 CONTINUE 

341 CONTINUE 

DO 370 1=1,3 

R4F101 I ) =TSG240*RD10F { I ) 

R4F 1 2 ( I ) =TSC240*RD12F( I ) 

R4F13 ( I ) =TSG240*RD13F ( I ) 

RD4F10I I )=TAU45*RCD10F( I ) 

RD4F12I I )=TAU45*KCD12F( I ) 

RD4F 1 3 ( I ) = TAU45»RCD13F( I ) 

370 CONTINUE 

CX10 = DSQRT ( (R4F10 ( 1 ) *»2 ) + ( R4F10 < 2 )*«2 ) + < R4F10 ( 3) «»2 ) ) 

380 DO 390 1=1,3 

X4F10( I )=XAST10( I J+R4F10 t I) 

X4F 1 2 ( I)=XAST12( I ) +K4F 12 ( I ) 

X4F131 I)=XAST13( I 1+R4F131 I ) 

X4F10I I + 3)=XAST10( I + 3)+RC4F10I I ) 

X4F 12 ( I+3)=XAST12( I+3)+RC4F12( I ) 

X4F13I I+3)=XAST13( I+3)+RC4F13( I ) 

390 CONTINUE 
69 DO 18 J= 1 , 3 

Y10(J)=X4F10( J)*CML 
Y12tJ)=X4F12( J)»CML 
Y 1 3 I J ) = X4F 1 3 ( J ) *C ^ L 
Y10(J+3)=X4F10(J+3)*Cl'W 
Y 1 2 ( J+3) =X4F12( J+3MCNV 
Y13(J+3)=X4F13(J+3)»C^V 
18 CONTINUE 

R10 C 2=Y1C(1)**2+Y10(2)*«2+Y1C(3)«*2 

R20 c 2= ( Y1C( 1 )-Y12(l) ) **2+ ( Y1C12 ) -Y12I 2) )**2+( Y10(3)-Y12< 3) )*«2 



R10=CSCRT(R10E2) 

R20=DSCRT I R20E2 ) 

IF! TPR16027, 6026, 6027 
6027 CONTINUE 

TAUP=4.0CC0«TAU 
WRITE ( 3, 6CC0) R20 

6000 FOR w AT ( I 1H R20IKM) = ,D16.8) 

WR I TE ( 3 » 6C01 ) TAUP 

6001 FOR M AT ( 11H FOUR TAU =,016.8) 

WRITEI3.6C02) CX10 

6002 FORMAT ( 24H MAG GF RIO CORRECTION =,C16.8> 

WRITE ( 3 , 6C05) DELT 

6005 FOR M AT ( 6E DEL T= , 0 1 6 . 0 > 

6026 CONTINUE 

I F ( SWT— 1.CC0J6024, 1C 01, 6 02 A 
6024 R10F4T=DABS(DELT) 

IF (R10F4T 16021, 6023,6021 

6021 R10F4T=AT/R10F4T 

I F ( 0 10F4T - 6.000)6022,6022,6023 

6022 TAU=TAU*.5CC0* (R10F4T+1 .000) 

1006 IFIONEDAY- ( T AUN+ A » T AU* 1 . 1 DCO ) ) 1C CO, 1000,9 

1000 TAU=(ONEL)AY-TAUN) / A 
SWT=1.D0C 

GO TO 9 

1001 TOTD A Y=TCTDA Y+ONE D AY 
TDAYCU=(TCTDAY«CMT)/864CC.CCC 
SWT=O.DOG 

TIMF0=T0AYCU 

WRITEI3, 2C ) (Y10( I ), 1 = 1,6) .TIMED 
WRITE (3, 19 ) ( Y12( I ) , 1 = 1 , 6 ) 

WRITE (3, 19) ( Y13I I ) , 1=1,6) 

19 FORMAT ( IE 1 ,.6016.8) 

20 FORMAT ( 1E0, 7016. 8 ) 

2050 FORMA T ( 1 H 3D24.16) 

WRITE (3, 28)CX10,R10,R20,TIMEC,C0LNT 
28 FORMATQH ,5016.8) 

IF( TIMEMX-TDAYCU) 1C03, 1C03, 1C02 

1002 CONTINUE 
TAUN=O.OCO 
GO TO 6024 

6023 R10F4T=6.CC0 
GO TO 6022 

1003 i f i more . eg . i ) go to 900 

65 CONTINUE 
STO D 
END 

$ I B FTC SU R 1 LIST, REF, NCDECK 
SUBROUTINE SUB1 (X,M) 

DOUPLE PRECISION 
0 X ( 6 ) , M , 

DERR.TIMEDN, ITDN, IT0,A,CX10, 

DX10.X12, X13,XTEMP,XKEP,XCEL, 

DST2X10, ST2X12 , S T2 X 1 3 , S T 2 X20 , ST2X 30, ST 2X 23 , 
DM1,"2,M3,M12,M13,M23,M2F12,M2F23,M3F13,M3F23, 

DCML,CMV,CMT, Y10.Y12, Y13.TIMEC.TI MEMX, 

DTAUC, TAU,R10,R20,TEMP,Q,R1C£2,R2CE2, 

OL, Cl, C2,C3, DELTA, H,XFR1G,XAFR10, 

OXFR12,XFR13,XFR20,XER23,XFR30,XAFR12,XAFR13,XAFR20,XAFR23,XAFR30» 
DRE2,R,KS I,ETA,ZETA,ChI , R 1 2 , R 1 3 , R 2 3 , R 30, Z , 

DR10F2C,R2CE2C,R10C,R20C 



DOUBLE PRECISION fc V. 1C , W W 1 2 , Wfe 13 , 1- fc 2 C , WW30 , V, W2 3 , 
DFF1C,FF12,FF13,FF2C,FF3C,FF23 
OOUPLE PRECISION 

OW.TCTDAY ,TDAYCU,ONEDAY, SWT.TAUN, COUNT, AT, BT, 

DRD10,RC1 2,RD13,RDC10,RDD12,RCD13 , 

0DEL10»DEL12»DEL13,MCNE*X20»0EL20 »X 30 , DE L30 , X2 3 , 0EL23 , 
0EPS10,XAST10,EPS12,XAST12,EPS13,XAST13,XAST20,XAST23,XAST30, 
DRAST10,RAST12,RAST13,RAST20,RAST23,RAST30,S10,S12,S13,S20, 
0S23.S30* R00T10,R00T12,RDCT13,RD1CF, RD1 2F , RD 1 3F , RC010F , RDD 1 2F , 
DKCD13F,RWGT,RDWGT , TSQ240 , R4F 10 , R4F l 2 , R4F 13 , TAU4 5 ,RD4F 10 , 
08X1C,BX12,BX13, 

DRD4B12*RD4F13»X4F10»X4FI2»X4F13 
COM VON XIC(12),X12(12),X13(12),XTEMP(6) , XK EP ( 6 ) , XCEL ( 6 ) , 
CM1,«2,M3,M12,M13,M23, M2F12, M2F23, M3F13.V3F23, 
CST2X1C(6>,ST2X12(6),ST2X13(6),ST2X20(6) , ST2X3C < 6 ) , ST2X23 ( 6) , 
CCML ,C«V,CMT,Y1Q(6),Y12(£),Y13(6) .TIMED, TIMEMX, 

CTAUC, TAU.R10, R20 , TEMP , Q , R 1 C 62 , R2 CE 2 , 

CL,C1,C2,C3,CELTA,H,XFR1C ( 3 ) , X AF R 10 ( 3 ) , 

CRE2,R,KSI ,ETA,ZETA,CHI ,R12< 6 ),R13(6),R23(6),R30(6),Z, 
CR1QF2C,R20E2C,R10C,K20C 

COM VON WW10(6) ,V.Wl2( 6) ,V*W13(6) ,Wk20(b) , V>W30<6) ,WW23(6) , 

CFF1CI6) , FF 12 l 6 ) ,FF13(6) ,FF20(6) ,FF3C(6) » FF2 3 ( 6 ) 

COM VUN 

CW, TOTDAY, TCAYCU.ONEDAY, SWT , TAUN , COUNT , AT,BT, 

CTSQ 240,8X10 (6) ,8X12(6), 8X13(6), 

CRD012FI6 ) ,RC013F( 6) ,RV.GT (4) , RDWG T < 4 ) , R4 F 10 1 3 ) , R4F 1 2 ( 3 ) , 

CRD10I3), RD 12(3), RD 1313), RCC10(3),RCD12(3), ROD 13(3), 

CQEL10 ( 6 ) , CEL 1 2 ( 6 ) «DtL13(6) ,DEL20(6) ,DEL23(6) ,DEL30(6), 

CMON p , X 20 ! 6 ) ,X23(b),X3C(6),EPS10(6),£PS12(6),tPS13(6), 
CXAST10(6),XAST12(6),XAST13(6),XAST2C(6),XAST23(6),XAST30(6), 
CRASTlOIb ) .RAST121 6) ,RAST13(b),RAST20(6) ,RAST23 ( 6) ,RAST30( 6) , 
CS10(6)»S12(6)»S13(6)»S20(6),S23(6),S30(fc) , R DOT 10(6) , 

CRD0T12(6 ) ,RC0T13( 6) ,KClCF(6),RD12F(6),RD13F(b),RCD10F(6), 

CR4F 13(3) ,TAU45,KD4F10(3) ,RC4F12(3) , RD4F 13 ( 3 ) , X4F10 ( 6 ) ,X4F12<6), 

C X4F 1 3 < 6 ) »XFR1 2 ( 3 ) ,XFR13(3) ,XFR20 (3),XFR23(3),XFR30(3),XAFR12(3), 
C X AF R 1 3 ( 3 ) , X AFR20 ( 3 ) , XAFR23(3),XAFR3C(3) 

RE2=X( l) *«2+X(2)*«2+X(3)**2 
R=DSQRT ( RE2 ) 

KSI=M*( T AU**2 ) /R**3 

ETA=(X(l)»X(4)+X(2)*X(5)+X(3)*X(fc) )*TAU/RE2 

ZETA=(X(4)**2+X(3)**2+X(fc)**2)»(TAU**2)/RE2-KSI 

CHI=KSI-ZETA 

Z=MCNE-.5CC0*ETA-( l.DCOZ 6 . DOO ) * ZE T A+ . 5C00*E T A « * 2+ ( 5 . 000/ 1 2 . 000 ) 
1*ETA*ZETA 
I TER=0 
I TER= I TER+ 1 
L = CH I * 1**2 

C2=.5D00*(MGN£-(L/12.D00)*(MONE-(L/30.DC0)* 

1 (M0N£-(L/56.000)*(M0NE-(L/90.CC0> ) ) ) ) 

C 3= ( MONE / b.DOO)*(MONE-(L/2G.OOO )< ( MONE- (L/42.DC0)* 

1 ( MONE- (L/72.DC0)* (MONE- (L/11C.DC0) ) ) ) ) 

C1 = M 0NE-L*C 3 

DELTA=M0NE+C1*ETA*Z+C2*ZETA*Z**2 
H = C?*£TA*Z«*2 + C3»ZETA*Z**3 + Z-M,0N£ 

Z=Z-H/CELTA 

IF(DABS(F).LE. 1.0-07) GO TC 31 
IF( ITER. LE. 10) GO TO 30 
CONTINUE 
DO 40 1=1,3 

XDEL ( I)=(-KSI*Z»*2)*IC2*X(I ) +C3« Z *TAU*X < I +3 ) ) 



XDEL(I+3)=(-KSI*£)*(Cl*X(I)+C2*Z«TAU*X( 1+3) )/(DELTA*TAU) 
XKE P ( I ) =X ( I ) + T AU* X ( 1 + 3 ) + XD E L ( I ) 

XKEPU+3) = XC I + 3 ) + XDE L { 1+3) 

40 CONTINUE 
RETURN 
END 

$DATA 


• 2 3 62 1 5 5 1 D06 
.84415966DCC 
•22841895D06 
.8243782IDC0 
-.3I9572I6D08 
-.2 862*154002 
7 • 5 D- 08 


284078150C6 
-.715041670-1 
-.286I2319DC6 
. 50846563CC0 
. L 3 642 L31D09 
56407282DC1 
1.D-C2 


- . 1 594 4 7 98 CO 6 
- . 4 C 7 8 1 9 76 C- 1 
-• 1 60 7 82 32 CQ6 
• 18637202 COO 
.59157903008 
-.2*463687001 
1.0C01 


17 9. DOO 
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